Time varying resistance coronary boundary condition

docs/references.bib

@@ -126,6 +126,18 @@ @article{caruel13
  year = {2013}
 }
 
+@article{yong25,
+  title = {A lumped parameter model of the coronary circulation incorporating time-varying resistance, intramyocardial pressure and vascular compliance},
+  author = {Enhui Yong and Javeria Latief and Yufei Wang and David Erlinge and Axel Dahlgren and Tushar Kotecha and Vivek Muthurangu and Ryo Torii},
+  journal = {Journal of Biomechanics},
+  volume = {189},
+  pages = {112679},
+  year = {2025},
+  doi = {10.1016/j.jbiomech.2025.112679},
+  url = {https://doi.org/10.1016/j.jbiomech.2025.112679},
+  publisher = {Elsevier}
+}
+
 @article{Regazzoni2022,
 title = {A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation},
 journal = {Journal of Computational Physics},
@@ -139,4 +151,3 @@ @article{Regazzoni2022
 keywords = {Mathematical modeling, Multiphysics models, Multiscale models, Cardiac modeling, Cardiac electromechanics},
 abstract = {We propose a novel mathematical and numerical model for cardiac electromechanics, wherein biophysically detailed core models describe the different physical processes concurring to the cardiac function. The core models, once suitably approximated, are coupled by a computationally efficient strategy. Our model is based on: (1) the combination of implicit-explicit (IMEX) schemes to solve the different core cardiac models, (2) an Artificial Neural Network based model, that surrogates a biophysically detailed but computationally demanding microscale model of active force generation and (3) appropriate partitioned schemes to couple the different models in this multiphysics setting. We employ a flexible and scalable intergrid transfer operator, which allows to interpolate Finite Element functions between nested meshes and, possibly, among arbitrary Finite Element spaces for the different core models. Our core 3D electromechanical model of the left ventricle is coupled with a closed-loop 0D model of the vascular network (and the other cardiac chambers) by an approach that is energy preserving. More precisely, we derive a balance law for the mechanical energy of the whole circulatory network. This provides a quantitative insight into the energy utilization, dissipation and transfer among the different compartments of the cardiovascular network and during different stages of the heartbeat. On this ground, a new tool is proposed to validate some energy indicators adopted in the daily clinical practice. A further contribution of this paper is the proposition of a robust algorithm for the reconstruction of the stress-free reference configuration. This feature is fundamental to correctly initialize our electromechanical simulations. As a matter of fact, the geometry acquired from medical imaging typically refers to a configuration affected by residual internal stresses, whereas the elastodynamics equations that govern the mechanics core model are related to a stress-free configuration. To prove the biophysical accuracy of our computational model, we address different scenarios of clinical interest, namely by varying preload, afterload and contractility.}
 }
-

scripts/OpenLoopCoronaryBC.yaml

@@ -0,0 +1,27 @@
+variables:
+  - Pin
+  - Qin
+  - Vim
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dVim_dt
+
+constants:
+  - Ra
+  - Ram
+  - Rv
+  - Ca
+  - Cim
+  - Pim
+  - Pv
+  - Pim_0
+  - P_Cim_0
+
+time_dependent:
+  - Pim
+
+residuals:
+  - Cim * Rv * Qin - Vim + Cim * (-P_Cim_0 + Pim_0 - Pim + Pv) - Cim * Rv * dVim_dt - Ca * Cim * Rv * dPin_dt + Ra * Ca * Cim * Rv * dQin_dt
+  - Cim * Rv * Pin - Cim * Rv * Ra * Qin - Rv * Vim - Cim * (Rv + Ram) * (P_Cim_0 - Pim_0 + Pim) - Cim * Rv * Ram * dVim_dt - Ram * Vim + Ram * Cim * Pv

scripts/OpenLoopCoronaryVarResBC.yaml

@@ -0,0 +1,28 @@
+variables:
+  - Pin
+  - Qin
+  - Vim
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dVim_dt
+
+constants:
+  - Ra
+  - Ram
+  - Rv
+  - Ca
+  - Cim
+  - Pim
+  - Pv
+  - Pim_0
+  - P_Cim_0
+
+time_dependent:
+  - Pim
+  - Ram
+
+residuals:
+  - Cim * Rv * Qin - Vim + Cim * (-P_Cim_0 + Pim_0 - Pim + Pv) - Cim * Rv * dVim_dt - Ca * Cim * Rv * dPin_dt + Ra * Ca * Cim * Rv * dQin_dt
+  - Cim * Rv * Pin - Cim * Rv * Ra * Qin - Rv * Vim - Cim * (Rv + Ram) * (P_Cim_0 - Pim_0 + Pim) - Cim * Rv * Ram * dVim_dt - Ram * Vim + Ram * Cim * Pv

scripts/jacobian.py

@@ -25,13 +25,11 @@ def load_model(filepath):
     time_dependent = {context[s] for s in data.get('time_dependent', [])}
 
     assert len(variables) == len(derivatives), f"Number of variables must be equal to number of derivatives"
-    assert len(variables) - 2 == len(residuals), f"Number of residuals must be number of unknowns minus 2"
 
     return variables, derivatives, constants, residuals, time_dependent
 
 def extract_linear(residuals, y, dy):
-    def is_constant_coeff(coeff):
-        return (coeff and coeff.free_symbols.isdisjoint(y.free_symbols | dy.free_symbols))
+    zero_subs = [(v, 0) for v in list(y) + list(dy)]
     nr = residuals.shape[0]
     ny = y.shape[0]
     E = Matrix.zeros(nr, ny)
@@ -40,8 +38,9 @@ def is_constant_coeff(coeff):
         for j in range(ny):
             for mat, dd in zip([E, F], [dy, y]):
                 coeff = residuals[i].coeff(dd[j])
-                if is_constant_coeff(coeff):
-                    mat[i, j] = coeff
+                const_coeff = simplify(coeff.subs(zero_subs))
+                if const_coeff:
+                    mat[i, j] = const_coeff
     return E, F
 
 def extract_nonlinear(residuals, E, F, y, dy):
@@ -53,19 +52,27 @@ def extract_nonlinear(residuals, E, F, y, dy):
 def depends_on(expr, symbols_set):
     return any(sym in expr.free_symbols for sym in symbols_set)
 
-def partition_terms(E, F, C, dC_dy, dC_dydot, time_dependent_symbols):
+def partition_terms(E, F, C, dC_dy, dC_dydot, variables, time_dependent_symbols):
     parts = {"constant": [], "time": [], "solution": []}
-    for i in range(E.shape[0]):
-        parts["solution"].append(("C", i, C[i]))
-        for j in range(E.shape[1]):
-            for label, mat in [("E", E), ("F", F), ("dC_dy", dC_dy), ("dC_dydot", dC_dydot)]:
-                if label.startswith('d'):
-                    target = "solution"
-                elif depends_on(mat[i, j], time_dependent_symbols):
-                    target = "time"
+    variable_symbols = variables.free_symbols
+    for label, mat in [("E", E), ("F", F), ("C", C), ("dC_dy", dC_dy), ("dC_dydot", dC_dydot)]:
+        for i in range(mat.rows):
+            for j in range(mat.cols):
+                if mat.cols == 1:
+                    expr = mat[i]
                 else:
-                    target = "constant"
-                parts[target].append((label, i, j, mat[i, j]))
+                    expr = mat[i, j]
+                if expr != 0:
+                    if depends_on(expr, variable_symbols):
+                        target = "solution"
+                    elif depends_on(expr, time_dependent_symbols):
+                        target = "time"
+                    else:
+                        target = "constant"
+                    if mat.cols == 1:
+                        parts[target].append((label, i, expr))
+                    else:
+                        parts[target].append((label, i, j, expr))
     return parts
 
 def get_expressions(parts, type):
@@ -125,7 +132,7 @@ def main(yaml_path):
     C, dC_dy, dC_dydot = extract_nonlinear(residuals, E, F, y, dy)
 
     # split into constant, time dependent and solution parts
-    parts = partition_terms(E, F, C, dC_dy, dC_dydot, time_dependent)
+    parts = partition_terms(E, F, C, dC_dy, dC_dydot, Matrix(list(y) + list(dy)), time_dependent)
 
     # print C++ code
     for section in ['constant', 'time', 'solution']:

src/model/BlockType.h

@@ -31,7 +31,8 @@ enum class BlockType {
   chamber_sphere = 15,
   blood_vessel_CRL = 16,
   piecewise_cosine_chamber = 17,
-  piecewise_valve = 18
+  piecewise_valve = 18,
+  open_loop_coronary_var_res_bc = 19
 };
 
 /**

src/model/CMakeLists.txt

@@ -23,6 +23,7 @@ set(CXXSRCS
   Model.cpp
   Node.cpp
   OpenLoopCoronaryBC.cpp
+  OpenLoopCoronaryVarResBC.cpp
   Parameter.cpp
   PressureReferenceBC.cpp
   ResistanceBC.cpp
@@ -52,6 +53,7 @@ set(HDRS
   Model.h
   Node.h
   OpenLoopCoronaryBC.h
+  OpenLoopCoronaryVarResBC.h
   Parameter.h
   PressureReferenceBC.h
   ResistanceBC.h

src/model/Model.cpp

@@ -21,6 +21,7 @@ Model::Model() {
        block_factory<ClosedLoopHeartPulmonary>()},
       {"ClosedLoopRCR", block_factory<ClosedLoopRCRBC>()},
       {"CORONARY", block_factory<OpenLoopCoronaryBC>()},
+      {"CORONARY_VAR_RES", block_factory<OpenLoopCoronaryVarResBC>()},
       {"FLOW", block_factory<FlowReferenceBC>()},
       {"NORMAL_JUNCTION", block_factory<Junction>()},
       {"PRESSURE", block_factory<PressureReferenceBC>()},

src/model/Model.h

@@ -31,6 +31,7 @@
 #include "Junction.h"
 #include "Node.h"
 #include "OpenLoopCoronaryBC.h"
+#include "OpenLoopCoronaryVarResBC.h"
 #include "Parameter.h"
 #include "PiecewiseCosineChamber.h"
 #include "PiecewiseValve.h"

src/model/OpenLoopCoronaryBC.cpp

@@ -2,17 +2,19 @@
 // University of California, and others. SPDX-License-Identifier: BSD-3-Clause
 #include "OpenLoopCoronaryBC.h"
 
+#include "Model.h"
+
 void OpenLoopCoronaryBC::setup_dofs(DOFHandler& dofhandler) {
   Block::setup_dofs_(dofhandler, 2, {"volume_im"});
 }
 
 void OpenLoopCoronaryBC::update_constant(SparseSystem& system,
                                          std::vector<double>& parameters) {
   auto Ra = parameters[global_param_ids[0]];
-  auto Ram = parameters[global_param_ids[1]];
-  auto Rv = parameters[global_param_ids[2]];
-  auto Ca = parameters[global_param_ids[3]];
-  auto Cim = parameters[global_param_ids[4]];
+  auto Ram = get_Ram(parameters, 0.0);
+  auto Rv = parameters[global_param_ids[1]];
+  auto Ca = parameters[global_param_ids[2]];
+  auto Cim = parameters[global_param_ids[3]];
 
   if (steady) {
     // Different equations for steady initial condition
@@ -27,23 +29,27 @@ void OpenLoopCoronaryBC::update_constant(SparseSystem& system,
     system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -1.0;
     system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = Cim * Rv;
     system.F.coeffRef(global_eqn_ids[1], global_var_ids[1]) = -Cim * Rv * Ra;
-    system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -(Rv + Ram);
 
     system.E.coeffRef(global_eqn_ids[0], global_var_ids[0]) = -Ca * Cim * Rv;
     system.E.coeffRef(global_eqn_ids[0], global_var_ids[1]) =
         Ra * Ca * Cim * Rv;
     system.E.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -Cim * Rv;
-    system.E.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -Cim * Rv * Ram;
   }
 }
 
+double OpenLoopCoronaryBC::get_Ram(std::vector<double>& parameters,
+                                   double time) const {
+  // Base class returns constant resistance (time parameter unused)
+  return parameters[global_param_ids[6]];
+}
+
 void OpenLoopCoronaryBC::update_time(SparseSystem& system,
                                      std::vector<double>& parameters) {
-  auto Ram = parameters[global_param_ids[1]];
-  auto Rv = parameters[global_param_ids[2]];
-  auto Cim = parameters[global_param_ids[4]];
-  auto Pim = parameters[global_param_ids[5]];
-  auto Pv = parameters[global_param_ids[6]];
+  auto Ram = get_Ram(parameters, model->time);
+  auto Rv = parameters[global_param_ids[1]];
+  auto Cim = parameters[global_param_ids[3]];
+  auto Pim = parameters[global_param_ids[4]];
+  auto Pv = parameters[global_param_ids[5]];
 
   if (steady) {
     system.C(global_eqn_ids[1]) = Pv;
@@ -53,6 +59,8 @@ void OpenLoopCoronaryBC::update_time(SparseSystem& system,
     system.C(global_eqn_ids[1]) =
         (Ram * Cim * Pv) -
         Cim * (Rv + Ram) * (Pim + this->P_Cim_0 - this->Pim_0);
+    system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -(Rv + Ram);
+    system.E.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -Cim * Rv * Ram;
   }
 }
 
@@ -63,8 +71,8 @@ void OpenLoopCoronaryBC::setup_initial_state_dependent_params(
   auto P_in_dot = initial_state.ydot[global_var_ids[0]];
   auto Q_in_dot = initial_state.ydot[global_var_ids[1]];
   auto Ra = parameters[global_param_ids[0]];
-  auto Ram = parameters[global_param_ids[1]];
-  auto Ca = parameters[global_param_ids[3]];
+  auto Ram = get_Ram(parameters, 0.0);
+  auto Ca = parameters[global_param_ids[2]];
   // Pressure proximal to Ca and distal to Ra
   auto P_Ca = P_in - Ra * Q_in;
   auto P_Ca_dot = P_in_dot - Ra * Q_in_dot;
@@ -73,5 +81,5 @@ void OpenLoopCoronaryBC::setup_initial_state_dependent_params(
   // Pressure proximal to Cim/Vim and distal to Ram
   this->P_Cim_0 = P_Ca - Ram * Q_am;
   // Initial intramyocardial pressure
-  this->Pim_0 = parameters[global_param_ids[5]];
+  this->Pim_0 = parameters[global_param_ids[4]];
 }

src/model/OpenLoopCoronaryBC.h

@@ -73,12 +73,12 @@
  * Parameter sequence for constructing this block
  *
  * * `0` Ra: Small artery resistance
- * * `1` Ram: Microvascular resistance
- * * `2` Rv: Venous resistance
- * * `3` Ca: Small artery capacitance
- * * `4` Cim: Intramyocardial capacitance
- * * `5` Pim: Intramyocardial pressure
- * * `6` Pv: Venous pressure
+ * * `1` Rv: Venous resistance
+ * * `2` Ca: Small artery capacitance
+ * * `3` Cim: Intramyocardial capacitance
+ * * `4` Pim: Intramyocardial pressure
+ * * `5` Pv: Venous pressure
+ * * `6` Ram: Microvascular resistance
  *
  * ### Usage in json configuration file
  *
@@ -124,13 +124,13 @@ class OpenLoopCoronaryBC : public Block {
       : Block(id, model, BlockType::open_loop_coronary_bc,
               BlockClass::boundary_condition,
               {{"Ra1", InputParameter()},
-               {"Ra2", InputParameter()},
                {"Rv1", InputParameter()},
                {"Ca", InputParameter()},
                {"Cc", InputParameter()},
                {"t", InputParameter(false, true)},
                {"Pim", InputParameter(false, true)},
                {"P_v", InputParameter()},
+               {"Ra2", InputParameter()},
                {"closed_loop_outlet", InputParameter(true, false, false)}}) {}
 
   /**
@@ -179,9 +179,42 @@ class OpenLoopCoronaryBC : public Block {
    */
   TripletsContributions num_triplets{5, 4, 0};
 
+  /**
+   * @brief Virtual destructor for proper cleanup in derived classes
+   */
+  virtual ~OpenLoopCoronaryBC() = default;
+
  protected:
+  /**
+   * @brief Protected constructor for derived classes
+   *
+   * Allows derived classes to specify their own block type and parameters
+   *
+   * @param id Global ID of the block
+   * @param model The model to which the block belongs
+   * @param block_type Type of the block
+   * @param input_params Input parameters for the block
+   */
+  OpenLoopCoronaryBC(
+      int id, Model* model, BlockType block_type,
+      const std::vector<std::pair<std::string, InputParameter>>& input_params)
+      : Block(id, model, block_type, BlockClass::boundary_condition,
+              input_params) {}
+
   double P_Cim_0 = 0;  ///< Pressure proximal to Cim/Vim at initial state
   double Pim_0 = 0;    ///< Pim at initial state
+
+  /**
+   * @brief Get microvascular resistance value
+   *
+   * For the base class, returns constant Ram. Derived classes may
+   * override to provide time-varying resistance.
+   *
+   * @param parameters Parameters of the model
+   * @param time Current simulation time (unused in base class)
+   * @return Microvascular resistance Ram
+   */
+  virtual double get_Ram(std::vector<double>& parameters, double time) const;
 };
 
 #endif  // SVZERODSOLVER_MODEL_OPENLOOPCORONARYBC_HPP_

src/model/OpenLoopCoronaryVarResBC.cpp

@@ -0,0 +1,38 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+#include "OpenLoopCoronaryVarResBC.h"
+
+#include <cmath>
+
+#include "Model.h"
+
+double OpenLoopCoronaryVarResBC::get_Ram(std::vector<double>& parameters,
+                                         double time) const {
+  // Extract parameters for time-varying resistance
+  auto Ram_min = parameters[global_param_ids[6]];
+  auto Ram_max = parameters[global_param_ids[7]];
+  auto T_vc = parameters[global_param_ids[8]];
+  auto T_vr = parameters[global_param_ids[9]];
+
+  // Get time within current cardiac cycle
+  double t_cycle = fmod(time, model->cardiac_cycle_period);
+
+  // Compute e(t) based on phase in cardiac cycle
+  double e_t;
+  if (t_cycle <= T_vc) {
+    // Contraction phase
+    e_t = 0.5 * (1.0 - cos(M_PI * t_cycle / T_vc));
+  } else if (t_cycle <= T_vc + T_vr) {
+    // Relaxation phase
+    e_t = 0.5 * (1.0 + cos(M_PI * (t_cycle - T_vc) / T_vr));
+  } else {
+    // Rest phase
+    e_t = 0.0;
+  }
+
+  // Compute time-varying resistance
+  double sqrt_Ram_min = sqrt(Ram_min);
+  double sqrt_Ram_max = sqrt(Ram_max);
+  double Ram_t = sqrt_Ram_min + (sqrt_Ram_max - sqrt_Ram_min) * e_t;
+  return Ram_t * Ram_t;
+}