<figure><img src="image_1.png" /></figure>

HAL
open science

# PhysioBlocks: an Opensource Python Library for Simulating Block Diagrams of Dynamical Physiological Systems

Dominique Chapelle, Colin Drieu, François Kimmig

► To cite this version:

Dominique Chapelle, Colin Drieu, François Kimmig. PhysioBlocks: an Opensource Python Library for Simulating Block Diagrams of Dynamical Physiological Systems. 2025. hal-05436182

HAL Id: hal-05436182

https://hal.science/hal-05436182v1

Preprint submitted on 31 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# PhysioBlocks: an Opensource Python Library for Simulating
Block Diagrams of Dynamical Physiological Systems

Dominique Chapelle<sup>a,b</sup>, Colin Drieu<sup>a,b</sup>, François Kimmig<sup>a,b</sup>

<sup>a</sup> Inria

<sup>b</sup> CMAP, Ecole Polytechnique, CNRS, Institut Polytechnique de Paris

## Abstract

The PhysioBlocks Python library is designed to simulate the dynamics of physiological systems (in particular cardiovascular systems) represented by block diagrams, in order to provide built-in modularity. Accordingly, a system is represented by a network of modules (blocks) connected by nodes in which they share quantities (degrees of freedom) and exchange fluxes. The user can easily create a new network by combining existing blocks. At a more advanced level, new blocks can be defined. The library is distributed under the LGPL-3.0-only license, and the initial distribution focuses on providing building blocks associated with lumped-parameter models of the cardiovascular system.

**Keywords:** Block diagram; Dynamical system; Physiology; Cardiovascular system

## 1 Introduction

The modeling and numerical simulation of dynamical physiological systems is a vast and extremely active domain, driven by a wealth of biological and medical applications, in particular. In this context, the initial motivations of the PhysioBlocks library stem from two main realizations:

1. Specific applications typically require a limited number of novel components, to be inserted within a network of other previously elaborated and tested building bricks;

2. Such modeling networks often take the form of a set of interconnected bricks that share some physical quantities (e.g. pressures, concentrations, ...) at the connection locations – that we call nodes – and exchange various types of fluxes at the same locations.

Focusing on dynamical systems that can be represented by this type of network, the block diagram formalism – heavily used in electrical engineering and control theory, in particular – is very natural, and also very powerful as an underlying principle for devising a software library that has the built-in modularity and versatility required, and this is the major paradigm of PhysioBlocks. This paradigm makes the design of new systems based on pre-existing blocks a very simple task that requires no coding, i.e., only the description of the new network. Moreover, incorporating a new physiological component in the library is also simplified in the sense that it requires no other coding than that defining the specific block associated with this new component. The library is primarily designed with so-called “lumped parameter models” in mind, i.e. for networks in which individual building blocks typically represent dynamical systems associated with a few ordinary differential equations (ODEs). Nevertheless, it can also be used with some much “heavier” building blocks (e.g., a 3D finite element model of the heart)
- probably relying on compiled code encapsulated within a Python wrapper – for which the rest of the network could be thought of as a boundary condition.

As a primary example in our presentation, we will use the cardiovascular system, and the building blocks provided in the initial distribution of the library are focused on this type of application, indeed. In particular, we include in the library the simplified heart model – or more strictly speaking, cavity model – that has been developed at Inria for over a decade, see initial description in [3]. We emphasize that this model is a “0D” – namely, consisting of a set of ODEs – reduction of the 3D biomechanical heart model presented in [5], for which the reduction process solely relies on geometric and kinematical approximations. This means that the tissue mechanics – and the multiscale contraction law in particular, based on the Huxley’57 description of actin-myosin interaction – is carefully preserved as in 3D.

This implies that the 0D and 3D models – together with a 1D model well-adapted to uni-axial in-vitro experiments with elongated samples such as papillary muscles or isolated car-diomyocytes, see [3] – can be thought of as members of a single hierarchical family, so that, for instance, the 0D model can be used to pre-calibrate a 3D model at a much reduced computa-tional cost. Nevertheless, the 0D model itself has already been used and validated in a variety of actual medical applications, in cardiology to start with [13, 7, 6], and also in critical care [11]. Therefore, we believe that this model can be of wide interest in other applications as well, and that sharing it as part of an opensource software is a valuable contribution to the biomedical community concerned with cardiovascular modeling. Of course, other popular cardiac lumped parameter models can be implemented and / or contributed later, such as the varying elastance model based on [14], or the so-called “single fiber model” of [1] implemented in the CircAdapt<fnref n="1" /> cardiovascular simulator.

Contributions of new blocks are strongly encouraged, regardless of whether they directly or indirectly concern the cardiovascular system, and especially if they satisfy the essential requirements of numerical stability prevailing in the library. After validation by the core team, such contributions will be distributed with the future versions of the library, and under the same LGPL-3.0-only license (while preserving the copyright associated with their actual authors).

The rest of the paper is organized as follows. Section 2 presents the major principles of the PhysioBlocks library, in adequation with the above-stated motivations. Then, Section 3 illustrates the detailed use of the library in the case of a simple – yet realistic – cardiovascular system, including the above-mentioned physiological cardiac cavity model. Section 4 gives some conclusions, together with important development directions. Finally, the appendices provide the list of all building blocks provided in the initial distribution (Appendix A), the detailed code for a simple $R-C$ Windkessel block (Appendix B), and eventually the detailed equations implemented in the cardiac cavity block (Appendix C).

# 2 Major Principles

## 2.1 Main concepts and objects

The type of systems that we consider in our block diagram formalism are represented by networks of elementary modules (blocks) connected by nodes at which they share quantities (degrees of freedom) and exchange fluxes. Accordingly, the main objects defined and manipulated in the library are the following:

* Block: individual modeling module

<footnote n="1">¹https://www.circadapt.org</footnote>
* **Node**: physical locations where **Blocks** are interconnected and share various physical quantities named **DOFs** (degrees of freedom), and where fluxes coming from connected modules satisfy balance equations (a.k.a. Kirchhoff laws), with one such equation per type of flux present at the node

* **Net**: overall specific system consisting of a set of **Blocks** connected via a set of **Nodes**

* **ModelComponent**: for some **Blocks** where the dynamics involves some variables other than the **DOFs** carried by the associated **Nodes**, the dynamics of these other internal variables can be optionally defined in such **ModelComponents**, for modularity purposes (in particular to share such equations between **Blocks** of different types)

In the example of the cardiovascular system, **Blocks** represent individual organs or compartments within the system, **DOFs** borne by the nodes are typically blood pressures, and the associated fluxes exchanged are blood flows.

## 2.2 Assembling and solution process

Based on the above-described architecture of **Nets**, the first stage of the solution process consists of looping over all the **Blocks** to obtain the equations governing all the unknowns at each time step, given the solution at the previous time step. In PhysioBlocks we favor implicit time schemes for stability purposes – in particular via the fulfillment of energy balances at the discrete level – and therefore with nonlinear physical laws – e.g. for a cardiac cavity – the corresponding governing equations are nonlinear. As a consequence, in the assembling process, the **Blocks** need to provide their contributions to the overall system of equations in terms of residuals and associated tangents to be used in a Newton-type algorithm. More specifically, each **Block**:

* must provide the equations defining all the relevant fluxes at every **Node** of the **Block** – i.e. the fluxes themselves and their tangents – which are all added into the equations representing the Kirchhoff laws for each flux at every **Node**, similarly to what is done in a finite element code to concatenate local stiffness operators into global ones, see Fig. 1;

* will directly provide the residuals and tangents for all equations governing its own internal variables, if any, possibly making use of **ModelComponents**.

<figure><img src="image_2.png" /><figcaption>Figure 1: System assembling example</figcaption></figure>

We emphasize that the fluxes defined in all the **Blocks** provided in the library should be understood as fluxes prevailing at time $t^{n+\frac{1}{2}} = \frac{t^n+t^{n+1}}{2}$ when seeking the solution for time $t^{n+1}$ via an implicit time scheme – again, for stability purposes via discrete-time energy balances – see example in the next subsection. Users implementing new **Blocks** should therefore aim at complying with this principle, in order to benefit from the numerical stability inherent to the pre-existing **Blocks**, and to avoid introducing consistency errors.
## 2.3 Example with simple Windkessel model

As a basic example, we consider a simple $R-C$ Windkessel model, see Figure in Section A.1.2. The pressures $P_1$ and $P_2$ are the **DOFs** at the **Nodes** of this **Block** – note that there are no internal variables in this **Block** – and in continuous time the fluxes are given by

$$Q_1 = \frac{P_2 - P_1}{R}, \quad Q_2 = \frac{P_1 - P_2}{R} - C\dot{P}_2. \qquad (1)$$

The corresponding mid-point expressions in discrete time are then

$$Q_1^{n+\frac{1}{2}} = \frac{P_2^{n+\frac{1}{2}} - P_1^{n+\frac{1}{2}}}{R}, \quad Q_2^{n+\frac{1}{2}} = \frac{P_1^{n+\frac{1}{2}} - P_2^{n+\frac{1}{2}}}{R} - C \frac{P_2^{n+1} - P_1^n}{\Delta t^n}, \qquad (2)$$

with $P_i^{n+\frac{1}{2}} = \frac{P_i^n + P_i^{n+1}}{2}$ for $i = 1, 2$, and $\Delta t^n = t^{n+1} - t^n$. This gives the contribution of the **Block** to the residuals associated with the Kirchhoff law for the fluxes at the **Nodes** and, since these expressions are linear with respect to the **DOFs**, the tangent expressions – i.e., the derivatives with respect to the **DOFs** – are straightforward.

## 2.4 User levels

There are three typical user levels for the library:

1. **Configure and run physiological systems simulations.** Here the user takes a pre-existing **Net**, and has the ability to adjust various parameters – e.g. to calibrate a model based on available data – before running simulations. This is a standard type of use in every simulation software, and the specificity of PhysioBlocks in this regard lies in the concept of *configuration aliases*, by which only the subset of parameters relevant to the user is being directly exposed, see Section 2.5.3 below for more detail.

2. **Create new systems with existing Blocks.** This is where the modularity paradigm of the library is paramount. Here, the user has the ability to build a new **Net without writing any code**, by connecting existing **Blocks** via a set of **Nodes**. The matching of **DOFs** and fluxes between all **Blocks** connected to a given **Node** is completely automatic based on their types as specified in the **Net** configuration, i.e. **DOFs** of the same type will coincide and fluxes of the same type will contribute to the same Kirchhoff law. This means that the user should carefully set these types when building a **Net** to ensure that the matching works as intended.

3. **Write and add new Blocks to the library.** Here, the user has the ability to devise a new **Block** when needed, based on the essential specifications prescribed in Section 2.2. Note that in some cases pre-existing **ModelComponents** may be reused in the new **Block**. As a simple example of **Block** construction, we give in Appendix B the code associated with the Windkessel model of Section 2.3.

## 2.5 Configuration files and aliases

Net layouts and simulation configurations are created via JSON files. While all configuration files share the same syntax, each user level will have a specific type of workflow with them.

1. **Configure simulations.** At Level 1, configuration files are used to update parameters for a given **Net**.
2. **Build Nets.** They allow a Level 2 user to create new **Net** layouts based on existing **Blocks** without writing new code.

3. **Write re-usable configuration subparts.** Through configuration aliases, a Level 2 or 3 user has a means to pre-configure technical parameters and write re-usable base configurations for higher user levels.

### 2.5.1 Net descriptions

The **Net** descriptions have two main components, see JSON Example 1:

1. **Nodes.** A list of the **Nodes** identifiers in the **Net**.

2. **Blocks.** A dictionary of all the **Blocks** in the **Net** associated with their identifiers. The **Blocks** are described by their:

(a) **Type.** The model type defining the flux and internal equations of the **Block**.

(b) **Flux type.** **Block** definitions being generic, the **Flux** type exchanged at its nodes is specified when creating a **Net**.

(c) **Global Ids.** The names in the **Net** of the **Block** parameters and variables (having local names in the **Block** definition); if global names are not specified, default names are used based on the **Block** Id and on the parameter or variable local name.

(d) **Connections.** For all the **Nodes** of the given **Block** type, their matching global **Node** Ids in the **Net**.

(e) **Submodels.** If the **Block** is enhanced with **ModelComponents**, they are described in this section (see JSON Example 3 for an illustration).

3. **DOF types definitions.** For each **Flux** type used in the **Net**, the name of the matching **DOF** type.

1 "nodes": [
2 "proximal",
3 "distal",
4 "venous"
5 ],
6 "blocks": {
7 "circulation_aorta_proximal": {
8 "type": "rc_block",
9 "flux_type": "blood_flow",
10 "nodes": {
11 "0": "proximal", # connection matching local node Id to global node Id
12 "1": "distal"
13 }
14 },
15 "circulation_aorta_distal": {
16 "type": "rc_block",
17 "flux_type": "blood_flow",
18 "nodes": {
19 "0": "distal",
20 "1": "venous"
21 }
22 }
23 },
24 "flux_dof_definitions": {
25 "blood_flow": "blood pressure"
26 }
JSON Example 1: Description of a simple circulation Net with two R-C Blocks

### 2.5.2 Simulation configuration

Simulation configuration files are used to launch a simulation instance. Their main components are (see JSON Example 2):

1. **Simulation type.** The type implementing the simulation process. Currently, we only provide a **Forward Simulation** type, that uses the solver (see next component) to find a solution at each time step, with an automatic strategy to adapt the time step duration if the solver does not converge. New simulation types will be implemented in future versions, like the **Estimation Simulation** using a Kalman filter for parameter estimation. The items used in the configuration file may evolve with the simulation type implementation.

2. **Solver.** The solver used for each simulation time step. Currently, we only provide a **Newton Solver** that implements a Newton method using the residual and tangents **System** described by the **Net**.

3. **Net.** The description of the **Net** to be used – most likely via an alias to an existing **Net** layout – into which we now incorporate adequate **Boundary Conditions** on some **Nodes**:

(a) **Flux** (Neumann boundary condition). The **Flux** is then added with other **Fluxes** of the same type shared at the **Node**.

(b) **DOF** (Dirichlet boundary condition). The **DOF** then becomes a parameter and the equation summing its matching **Flux** type is discarded from the global system.

4. **Time parameters.** The start time, duration and step size for the simulation.

5. **Variable initialization.** Their identifiers are determined from the **Net** description.

6. **Parameter values.** Their identifiers are also determined from the **Net** description. Their values can either be scalar, vector, boolean, string or function of time.
```json
{
  "type": "forward_simulation",
  "solver": {
    "type": "newton Solver"
  },
  "net": {
    "type": "circulation_alone_net",
    "boundary_conditions": {
      "proximal": [
        {
          "type": "condition",
          "condition_type": "blood_flow",
          "condition_id": "proximal.blood_flow"
        }
      ],
      "venous": [
        {
          "type": "condition",
          "condition_type": "blood pressure",
          "condition_id": "venous.blood pressure"
        }
      ],
      },
    },
    "time": {
      "start": 0.0,
      "duration": 25.0,
      "step": 1.0e-3
    },
    "variables_initializeation": {
      "proximal.blood pressure": 8.0e3,
      "distal.blood pressure": 8.0e3
    },
    "parameters": {
      "proximal.resistance": 8.0e6,
      "proximal.capacitance": 2.0e-10,
      "distal.resistance": 1.0e8,
      "distal.capacitance": 1.7e-08,
      "proximal.blood_flow": 0.5e-3,
      "venous.blood pressure": 1.6e3
    }
}
```

<caption>JSON Example 2: Configuration of a simulation for a simple circulation Net</caption>

In this example, we updated the **Net** description defined in JSON Example 1 with two boundary conditions. Note that the conditions at each **Node** are defined within a list because several **Flux** types can be shared at a single **Node**. We set a condition on the blood flow at the proximal **Node** and a condition on the blood pressure at the venous **Node**. Here, both parameters are set to constant values, but in a real use case we would typically use a time function to update the proximal blood flow for every time step.

### 2.5.3 Configuration aliases

Configuration aliases are designed to simplify JSON configuration files by defining re-usable configurations at various levels. Any configuration item can be defined in its own separate JSON file and will be reusable in any configuration. To use an aliased configuration, we call the alias file name or its full path in the type field of item configuration. Aliases are useful to:

1. **Define re-usable configuration items.** The JSON Examples 3 and 4 show how to alias a **Block** description. We define the alias in JSON Example 3 and name the file "spherical_cavity_hht.json" in a dedicated alias folder. We then use the aliased **Block** by
invoking its file name in a **Net** and complete its definition with its **Node** connections, see JSON Example 4.

```json
{
  "type": "cavity",
  "model_components": {
    "dynamics": {
      "type": "spherical_dynamics",
      ...
    },
    "velocity_law": {
      "type": "velocity_law_hht",
      ...
    },
    "rheology": {
      "type": "rheology_fiber_additive",
      ...
    },
    "model_components": {
      "active_law": {
        "type": "active_lawamiro_huxley_two_moments",
        ...
      }
    }
  }
}
```

<caption>JSON Example 3: Simplified **Block** alias definition for a spherical cavity</caption>

```json
{
  "cavity": {
    "type": "spherical_cavity_hht",
    "nodes": {
      "0": "ventricle"
    }
  }
}
```

<caption>JSON Example 4: Spherical cavity alias used in a **Net** description</caption>

2. **Pre-configure parameters.** In the JSON Example 5, we use the configuration from JSON Example 2 as an alias. Configuration items that use the same key as in the alias are overridden. Configuration items that are not overridden in the new definition keep the default value provided in the alias, and therefore they are hidden from a “high-level user” who only needs to update the configuration of JSON Example 5 to launch a simulation, making it much more user-friendly.

```json
{
  "type": "forward_simulation_circulation_alone",
  "time": {
    "start": 10.0,
    "duration": 60.0
  },
  "variables_initializeation": {
    "proximal.blood{j pressure": 8.0e3,
    "distal.blood pressure": 8.0e3
  },
  "parameters": {
    "proximal.resistance": 7.9e6,
    "distal.resistance": 1.2e8,
    "venous.blood pressure": 1595.0,
    "proximal.blood flow": 1.5e-3,
  }
}
```

<caption>JSON Example 5: Simplified configuration using an alias to pre-configure a simulation.</caption>
# 3 Example: Simple Cardiovascular System

## 3.1 Simple cardiovascular system net description

As an example, we build a **Net** configuration for the simple cardiovascular system depicted in Fig. 2. The associated **Block-Node** diagram is shown in Fig. 3.

<figure><img src="image_3.png" /><figcaption>Figure 2: Simple cardiovascular system with single cardiac chamber and arterial tree</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Simple cardiovascular system **Block-Node** representation</figcaption></figure>

The fluxes and internal equations implemented in each **Block** type and **ModelComponents** are detailed in Appendix A. An example of Python code for the *R-C Block* used in the **Net** is provided in Appendix B. Note that the **cavity Block** only contributes one flux depending on the cavity deformation. The **cavity Block** is enhanced with **ModelComponents** to implement – in a modular manner – the various equations governing the internal variables of the simplified cardiac model, see Appendix C.

We show the JSON description of the **Net** layout in JSON Example 6. Note that for the cardiac cavity, we used the alias previously defined in JSON Example 3.
1 "type": "net",
2 "nodes": [
3 "cavity",
4 "atrium",
5 "proximal",
6 "distal",
7 "venous"
8 ],
9 "blocks": {
10 "cavity": {
11 "type": "spherical_cavity_hht",
12 "flux_type": "blood_flow",
13 "nodes": {
14 "0": "cavity"
15 }
16 },
17 "valve_atrium": {
18 "type": "valve_r1_block",
19 "flux_type": "blood_flow",
20 "nodes": {
21 "0": "atrium",
22 "1": "cavity"
23 }
24 },
25 "valve_arterial": {
26 "type": "valve_r1_block",
27 "flux_type": "blood_flow",
28 "nodes": {
29 "0": "cavity",
30 "1": "proximal"
31 }
32 },
33 "capacitance_valve": {
34 "type": "c_block",
35 "flux_type": "blood_flow",
36 "nodes": {
37 "0": "cavity"
38 }
39 },
40 "circulation_aorta_proximal": {
41 "type": "rc_block",
42 "flux_type": "blood_flow",
43 "nodes": {
44 "0": "proximal",
45 "1": "distal"
46 }
47 },
48 "circulation_aorta_distal": {
49 "type": "rc_block",
50 "flux_type": "blood_flow",
51 "nodes": {
52 "0": "distal",
53 "1": "venous"
54 }
55 }
56 },
57 "flux_dof_definitions": {
58 "blood_flow": "blood Pressure"
59 }
## 3.2 Computational performances for the simple cardiovascular model

We measured the duration of 100 simulation runs for this model. The resulting system size is 11, and we set a 30 s simulation duration with time steps of 1 ms, recalling that we solve a nonlinear system by a Newton algorithm at each time step. We show the statistics of the simulation timings on a standard computer (M3 Max MacBook Pro) in the following table, where we see that we roughly achieve real time.

<table><thead><tr><td>Mean (s)</td><td>Median (s)</td><td>Minimum (s)</td><td>Maximum (s)</td></tr></thead><tbody><tr><td>29.77</td><td>29.82</td><td>27.98</td><td>31.14</td></tr></tbody></table>

## 4 Conclusions and Perspectives

The main motivation for this paper is to accompany the release of the initial version (V1) of the opensource PhysioBlocks library, of which we have described the major principles and ingredients.

Future distributions will incorporate new blocks, including for models other than that of the cardiovascular system, such as the respiratory system (with a view to cardiovascular-respiratory coupling).

Other important development directions concern:

* optimization of computational performances;

* incorporation of data assimilation features – satisfying the modularity principles advocated in [4] – in order to enable digital twin types of applications;

* development of a graphical user interface (GUI) to facilitate the design of new systems.

### Acknowledgements:

<figure><img src="image_5.png" /></figure>

The PhysioBlocks library is being developed as part of the MEDITWIN project (https://www.3ds.com/science/meditwin).

<figure><img src="image_6.png" /></figure>

This project has been funded by the French government as part of France 2030.

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

Funded by the European Union – Next Generation EU as part of the France Relance plan.

### References

[1] T. Arts, P.H. Bovendeerd, F.W. Prinzen, and R.S. Reneman. Relation between left ventricular cavity pressure and volume and systolic fiber stress and strain in the wall. *Biophysical Journal*, 59(1):93–102, 01 1991.

[2] J. Bestel, F. Clément, and M. Sorine. A biomechanical model of muscle contraction. In *MICCAI 2001, Lecture Notes in Computer Science*. Springer, 2001.
[3] M. Caruel, R. Chabiniok, P. Moireau, Y. Lecarpentier, and D. Chapelle. Dimensional reductions of a cardiac model for effective validation and calibration. *Biomechanics and Modeling in Mechanobiology*, 13(4):897-914, 12 2013.

[4] D. Chapelle, M. Fragu, V. Mallet, and P. Moireau. Fundamental principles of data assimilation underlying the Verdandi library: applications to biophysical model personalization within euHeart. *Medical & Biological Engineering & Computing*, 51(11):1221-1233, 11 2012.

[5] D. Chapelle, P. Le Tallec, P. Moireau, and M. Sorine. Energy-preserving muscle tissue model: formulation and compatible discretizations. *International Journal for Multiscale Computational Engineering*, 10(2), 2012.

[6] M. Gusseva, T. Hussain, C. Hancock Friesen, P. Moireau, A. Tandon, C. Patte, M. Genet, K. Hasbani, G. Greil, D. Chapelle, and R. Chabiniok. Biomechanical Modeling to Inform Pulmonary Valve Replacement in Tetralogy of Fallot Patients After Complete Repair. *Canadian Journal of Cardiology*, 37(11):1798-1807, 2021.

[7] M. Gusseva, T. Hussain, G. Hancock Friesen, C. and Greil, D. Chapelle, and R. Chabiniok. Prediction of Ventricular Mechanics After Pulmonary Valve Replacement in Tetralogy of Fallot by Biomechanical Modeling: A Step Towards Precision Healthcare. *Annals of Biomedical Engineering*, 49(12):3339-3348, 2021.

[8] P. Hauret and P. Le Tallec. Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact. *Computer Methods In Applied Mechanics And Engineering*, 195(37-40):4890-4916, 01 2006.

[9] H.M. Hilber, T.J.R. Hughes, and R.L. Taylor. Improved numerical dissipation for time integration algorithms in structural dynamics. *Earthq. Eng. Struct. Dyn.*, 5:283-292, 1977.

[10] F. Kimmig, D. Chapelle, and P. Moireau. Thermodynamic properties of muscle contraction models and associated discrete-time principles. *Advanced Modeling and Simulation in Engineering Sciences*, 6(1), 2019.

[11] A. Le Gall, F. Vallée, K. Pushparajah, T. Hussain, A. Mebazaa, D. Chapelle, E. Gayat, and R. Chabiniok. Monitoring of cardiovascular physiology augmented by a patient-specific biomechanical model during general anesthesia. A proof of concept study. *PLoS ONE*, 15(5), 2020.

[12] J. Manganotti, F. Caforio, F. Kimmig, P. Moireau, and S. Imperiale. Coupling reduced-order blood flow and cardiac models through energy-consistent strategies: modeling and discretization. *Advanced Modeling and Simulation in Engineering Sciences*, 8(1), 2021.

[13] B. Ruijsink, K. Zugaj, J. Wong, K. Pushparajah, T. Hussain, P. Moireau, R. Razavi, D. Chapelle, and R. Chabiniok. Dobutamine stress testing in patients with Fontan circulation augmented by biomechanical modeling. *PLoS ONE*, 15(2), 2020.

[14] H. Suga, K. Sagawa, and A.A. Shoukas. Load independence of the instantaneous pressure-volume ratio of the canine left ventricle and effects of epinephrine and heart rate on the ratio. *Circulation Research*, 32(3):314-322, 03 1973.
# A List of building blocks provided in initial distribution (V1)

## A.1 Blocks

### A.1.1 C Block

<figure><img src="image_9.png" /></figure>

$$ \text{Node 1:} \qquad Q_1 = -C\dot{P}_1 \qquad (3) $$

Discretised:

$$ Q_1^{n+\frac{1}{2}} = -C \frac{P_1^{n+1} - P_1^n}{\Delta t^n} \qquad (4) $$

### A.1.2 R-C Block

<figure><img src="image_10.png" /></figure>

$$ \text{Node 1:} \qquad Q_1 = \frac{P_2 - P_1}{R} \qquad (5) $$

$$ \text{Node 2:} \qquad Q_2 = \frac{P_1 - P_2}{R} - C\dot{P}_2 \qquad (6) $$

Discretised:

$$ Q_1^{n+\frac{1}{2}} = \frac{P_2^{n+\frac{1}{2}} - P_1^{n+\frac{1}{2}}}{R} \qquad (7) $$

$$ Q_2^{n+\frac{1}{2}} = \frac{P_1^{n+\frac{1}{2}} - P_2^{n+\frac{1}{2}}}{R} - C \frac{P_2^{n+1} - P_2^n}{\Delta t^n} \qquad (8) $$

### A.1.3 R-C-R Block

<figure><img src="image_11.png" /></figure>

$$ \text{Node 1:} \qquad Q_1 = \frac{P_{\text{mid}} - P_1}{R_1} \qquad (9) $$

$$ \text{Node 2:} \qquad Q_2 = \frac{P_{\text{mid}} - P_2}{R_2} \qquad (10) $$
**Internal equation:**

$$ \frac{P_1 - P_{mid}}{R_1} + \frac{P_2 - P_{mid}}{R_2} - C\dot{P}_{mid} = 0 \qquad (11) $$

**Discretized:**

$$ Q_1^{n+\frac{1}{2}} = \frac{P_{mid}^{n+\frac{1}{2}} - P_1^{n+\frac{1}{2}}}{R_1} \qquad (12) $$

$$ Q_2^{n+\frac{1}{2}} = \frac{P_{mid}^{n+\frac{1}{2}} - P_2^{n+\frac{1}{2}}}{R_2} \qquad (13) $$

$$ \frac{P_1^{n+\frac{1}{2}} - P_{mid}^{n+\frac{1}{2}}}{R_1} + \frac{P_2^{n+\frac{1}{2}} - P_{mid}^{n+\frac{1}{2}}}{R_2} - C \frac{P_{mid}^{n+1} - P_{mid}^{n}}{\Delta t^n} = 0 \qquad (14) $$

### A.1.4 R-L valve Block

<figure><img src="image_12.png" /></figure>

**Node 1:**

$$ Q_1 = -Q \qquad (15) $$

**Node 2:**

$$ Q_2 = Q \qquad (16) $$

**Internal equation:**

$$ L \dot{Q} + P_2 - P_1 + \begin{cases} RQ & \text{if } Q > 0 \\ R_{back}Q & \text{else} \end{cases} = 0 \qquad (17) $$

**Discretized:**

$$ Q_1^{n+\frac{1}{2}} = -Q^{n+\frac{1}{2}} \qquad (18) $$

$$ Q_2^{n+\frac{1}{2}} = Q^{n+\frac{1}{2}} \qquad (19) $$

$$ L \frac{Q^{n+1} - Q^n}{\Delta t^n} + P_2^{n+\frac{1}{2}} - P_1^{n+\frac{1}{2}} + \begin{cases} RQ^{n+\theta} & \text{if } Q^{n+\theta} > 0 \\ R_{back}Q^{n+\theta} & \text{else} \end{cases} = 0 \qquad (20) $$
### A.1.5 Cavity Block

<figure><img src="image_13.png" /></figure>

**Node 1:**

$$Q_1 = -\frac{dV(y)}{dt} \qquad (21)$$

with $V(y)$ given in (24).

**Discretization:**

$$Q_1^{n+\frac{1}{2}} = -\frac{V(y^{n+1}) - V(y^n)}{\Delta t^n} \qquad (22)$$

**Internal equations:** the dynamics of all internal variables – including the radial displacement $y$ – are provided by **Model-Components**, see below.

## A.2 Model components

The **ModelComponents** provided in the initial (V1) version of the library are those required to specify the dynamics of the internal variables of the above cardiac cavity **Block**. More specifically, we have four **ModelComponents** implementing the discrete equations given in Appendix C, namely,

* Velocity law: the **ModelComponent** named `velocity_law_hht` incorporates the two equations (29) and (30) into the model;

* Cavity dynamics: the **ModelComponent** named `spherical_dynamics` incorporates Eq. (31);

* Rheology: the **ModelComponent** named `rheology_fiber_additive` incorporates Eq. (32);

* Active law: the **ModelComponent** named `active_lawamiro_huxley_two_moments` incorporates the three equations Eqs. (33), (34) and (35).

## B Code for simple $R-C$ Windkessel block

```python
from dataclasses import dataclass
```

```python
import numpy as np
```

```python
import numpy as np
```

```python
from physioblocks.computing import Block, Expression, Quantity
```

```python
from physioblocks.computing import Block, Expression, Quantity
```

```python
from physioblocks.registers import register_type
```

```python
from physioblocks.simulation import Time
```

```python
# We register the class RCBlock with the id "rc_block" to be able
# to load and save a block instance in a configuration file.
```

```python
@register_type("rc_block")
```

```python
@dataclass
```

```python
class RCBlock(Block):
```

```python
# We are defining the blocks quantities here
```

```python
# We are defining the blocks quantities here
```

```python
# We are defining the blocks quantities here
```

```python
# We are defining the blocks quantities here
```
```python
pressure_1: Quantity np.float64 # Pressure at node 1 of the block
pressure_2: Quantity np.float64 # Pressure at node 2 of the block
resistance: Quantity np.float64 # Resistance value of the block
capacitance: Quantity np.float64 # Capacitor value of the block
time: Time # the simulation time

def flux_1(self) -> Any:
    """
   Computes the block flux out of node 1.
    """
    pressure_1 = 0.5 * (self.pressure_1.new + self.pressure_1.current)
    pressure_2 = 0.5 * (self.pressure_2.new + self.pressure_2.current)
    return (pressure_2 - pressure_1) / self.resistance.current


def dflux_1_dpressure_1(self) -> Any:
    """
   Computes the flux_1 partial derivative for pressure_1.
    """
    return - 0.5 / self.resistance.current


def dflux_1_dpressure_2(self) -> Any:
    """
   Computes the flux_1 partial derivative for pressure_2.
    """
    return 0.5 / self.resistance.current


def flux_2(self) -> Any:
    """
   Computes the block flux out of node 2.
    """
    pressure_1 = 0.5 * (self.pressure_1.new + self.pressure_1.current)
    pressure_2 = 0.5 * (self.pressure_2.new + self.pressure_2.current)
    delta pressure_2 = pressure_2.new - pressure_2.current
    return (
        (pressure_1 - pressure_2) / self.resistance.current
        - self.capacitance.current * self.time inv dt * delta
    )


def dflux_2_dpressure_1(self) -> Any:
    """
   Computes the flux_2 partial derivative for pressure_1.
    """
    return 0.5 / self.resistance.current


def dflux_2_dpressure_1(self) -> Any:
    """
   Computes the flux_2 partial derivative for pressure_2.
    """
    return (
        - 0.5 / self.resistance.current
        - self.capacitance.current * self.time inv dt
    )

# Define the expression for the flux going in the input node of the block
# with RCBlock methods
_rc_block_flux_1_expr = Expression(
    1, # the flux size
    RCBlock.flux_1, # the function to compute the flux
)
```
```R
76 "pressure_1": RCBlock.dflux_1_dpressure_1,
77 "pressure_2": RCBlock.dflux_1_dpressure_2,
78 }, # the functions to compute the flux partial derivatives
79 )
80
81 # Similarly define the expression for the flux going in the output node of
82 # the block
83 _rc_block_flux_2_expr = Expression(
84 1,
85 RCBlock.flux_2,
86 {
87 "pressure_1": RCBlock.dflux_2_dpressure_1,
88 "pressure_2": RCBlock.dflux_2_dpressure_2,
89 },
90 )
91
92 # link the flux expressions to the RCBlock Class
93 RCBlock.declares_flux_expression(
94 0, # Set the index carrying the flux
95 "pressure_1", # the local id of the dof associated with the flux
96 _rc_block_flux_1_expr # the flux expression
97 )
98 # Similarly for flux_2
99 RCBlock.declares_flux_expression(1, "pressure_2", _rc_block_flux_2_expr)
```

## C Detailed equations of cardiac cavity block

As originally proposed in [3], the geometry of the cardiac chamber – typically, the left ventricle – is approximated as a sphere, and all material properties are assumed to also satisfy spherical symmetry. This implies that the displacement field of the cardiac tissue enjoys spherical symmetry as well, and can be described by one scalar quantity $y$ representing the change of radius $R$ of the mid-surface of the wall with respect to the reference configuration, where the radius is denoted by $R_0$, namely, $R = R_0 + y$.

For the dynamical equation governing the displacement $y$, we consider the variant of the original equation of [3] proposed in [12] (see Eq. (10)-1 in the latter), i.e.

$$ \rho_0 \ddot{y} + \frac{\partial W_e(y)}{\partial y} + \frac{1}{R_0} \sum_v (y, \dot{y}) + \frac{k_s}{R_0} \left( \frac{y}{R_0} - e_c \right) = \frac{P_V}{|\Omega_0|} \frac{\partial V(y)}{\partial y}. \quad (23) $$

In this equation, the notation is defined as follows:

* $\rho_0$ denotes the tissue mass per unit volume in the reference configuration;

* $W_e(y)$ represents an elastic energy per unit volume (measured in the reference configuration) given by

$$ W_e(y) = C_0 \exp \left[ C_1 (2C(y) + C(y)^{-2} - 3)^2 \right] + C_2 \exp \left[ C_3 (C(y) - 1)^2 \right], $$

where $(C_0, C_1, C_2, C_3)$ denote material parameters of the passive constitutive law, and $C(y) = (1 + y/R_0)^2$ represents the component of the right Cauchy-Green deformation tensor along the fiber direction;

* $\Sigma_v(y, \dot{y})$ denotes a viscous stress defined by

$$ \Sigma_v(y, \dot{y}) = 2\eta (C(y) + 2C(y)^{-5}) \frac{\dot{y}}{R_0}, $$

with $\eta$ a solid viscosity parameter;
* $k_s$ denotes a passive elastic modulus of the sarcomeres, and $e_c$ the relative displacement of the myosin filaments with respect to the actin filaments (counted positively in the direction of sarcomere extension) divided by the length of a half-sarcomere;

* $P_V$ is the blood pressure inside the cavity;

* $|\Omega_0|$ is the total volume of cardiac tissue in the reference configuration, namely,

$$|\Omega_0| = \frac{4}{3}\pi \left[ \left(R_0 + \frac{d_0}{2}\right)^3 - \left(R_0 - \frac{d_0}{2}\right)^3 \right],$$

with $d_0$ the wall thickness in the reference configuration, and factorizing $R_0$ from the bracket a Taylor expansion readily gives the following approximation up to the third order in $d_0/2R_0$

$$|\Omega_0| = 4\pi R_0^2 d_0,$$

which is the expression that we will use in practice;

* $V(y)$ denotes the volume of the cardiac cavity, i.e.

$$V(y) = \frac{4}{3}\pi \left(R(y) - \frac{d(y)}{2}\right)^3 = \frac{4}{3}\pi \left(R_0 + y - \frac{d_0}{2}C(y)^{-1}\right)^3, \quad (24)$$

with $d(y)$ the wall thickness in the deformed configuration;

In addition, the equilibrium of forces within the sarcomere gives (see [10] for the rationale)

$$\mu \dot{e}_c + T_c = k_s \left( \frac{y}{R_0} - e_c \right), \quad (25)$$

where $\mu$ denotes a viscosity coefficient and $T_c$ the active tension resulting from the actin-myosin crossbridges.

Finally, the active tension $T_c$ is governed – together with an active stiffness quantity denoted by $K_c$ – by the following system of equations (given in [5], as an extension of the model proposed in [2] based on a moment approach applied with the Huxley'57 model),

$$\dot{K}_c = -(|u| + \alpha|\dot{e}_c|) K_c + n_0(e_c) K_0 |u|_+, \quad (26)$$

$$\dot{T}_c = -(|u| + \alpha|\dot{e}_c|) T_c + K_c \dot{e}_c + n_0(e_c) T_0 |u|_+, \quad (27)$$

where $u$ is a time-varying quantity representing electrical activation – positive during depolarisation / relaxation, and negative during repolarisation / contraction – $\alpha$, $K_0$ and $T_0$ are positive coefficients, and $n_0(e_c)$ is a positive function – with values typically lying between 0 and 1 – driving the fraction of maximum active tension and stiffness that can be reached in a sustained activation for a given sarcomere extension $e_c$, see [5]. Note that the typical behavior of $(K_c, T_c)$ resulting from this dynamical system is that it is driven towards $n_0(e_c)$ ($K_0, T_0$) when $u$ is being sustained at a positive level, indeed, and towards (0,0) when $u$ is being sustained at a negative value. In fact, for reasons pertaining to numerical stability we will also consider – as proposed in [5] – the auxiliary variable $\lambda_c = T_c/\sqrt{K_c}$, for which we infer from (26)-(27) the following dynamics

$$\dot{\lambda}_c = -\frac{1}{2}(|u| + \alpha|\dot{e}_c|) \lambda_c + \sqrt{K_c} \dot{e}_c + \frac{n_0(e_c)}{\sqrt{K_c}} \left( T_0 - \frac{K_0 \lambda_c}{2\sqrt{K_c}} \right) |u|_+. \quad (28)$$

We now present the time-discrete versions of the above equations as implemented in the library. For Eq. (23), we use the extended version of the so-called “HHT time integration scheme” of
[9] as proposed and analysed in [8], see also [12] for a specific analysis with the cavity model considered here, which establishes numerical stability with second-order accuracy. Accordingly, velocity and acceleration unknowns are introduced and solved for by

$$ \frac{\dot{y}^{n+1} - \dot{y}^n}{\Delta t^n} = (\frac{1}{2} + \alpha) \ddot{y}^{n+1} + (\frac{1}{2} - \alpha) \ddot{y}^n, \quad (29) $$

$$ \frac{y^{n+1} - y^n}{\Delta t^n} = \dot{y}^{n+\frac{1}{2}} + \frac{\alpha^2}{4} \Delta t^n (\ddot{y}^{n+1} - \ddot{y}^n), \quad (30) $$

with $\alpha$ some positive parameter. Eq. (23) is then discretized as follows

$$ \rho_0 \left[ \frac{\dot{y}^{n+1} - \dot{y}^n}{\Delta t^n} \right] + D W^{n+\frac{1}{2}\#} + \frac{1}{R_0} \Sigma_v (y^{n+\frac{1}{2}}, \dot{y}^{n+\frac{1}{2}}) + \frac{k_s}{R_0} \left( \frac{y^{n+\frac{1}{2}}}{R_0} - e_c^{n+\frac{1}{2}} \right) = \frac{P_V^{n+\frac{1}{2}}}{|\Omega_0|} D V^{n+\frac{1}{2}\#}, \quad (31) $$

with the notation

$$ D \Box^{n+\frac{1}{2}\#} = \begin{cases} \frac{\Box(y^{n+1}) - \Box(y^n)}{y^{n+1} - y^n}, & \text{if } y^{n+1} \neq y^n, \\ \frac{\partial \Box}{\partial y}(y^{n+\frac{1}{2}}), & \text{if } y^{n+1} = y^n = y^{n+\frac{1}{2}}. \end{cases} $$

The sarcomere equilibrium equation (25) is discretized using a standard mid-point rule

$$ \mu \frac{e_c^{n+1} - e_c^n}{\Delta t^n} - k_s \left( \frac{y^{n+\frac{1}{2}}}{R_0} - e_c^{n+\frac{1}{2}} \right) = -T_c^{n+\frac{1}{2}\#}, \quad (32) $$

with

$$ T_c^{n+\frac{1}{2}\#} = \lambda_c^{n+\frac{1}{2}} \sqrt{K_c^{n+1}}, \quad (33) $$

and for the discretization of ($K_c, \lambda_c$) we adopt the slight modification introduced in [12] of the original time scheme of [5] to use a mid-point – instead of backward-Euler – rule for $\lambda_c$, viz.

$$ \frac{K_c^{n+1} - K_c^n}{\Delta t^n} = - \left( |u^{n+1}| + \alpha \left| \frac{e_c^{n+1} - e_c^n}{\Delta t^n} \right| \right) K_c^{n+1} + n_0 (e_c^n) K_0 |u^{n+1}|_+, \quad (34) $$

$$ \frac{\lambda_c^{n+1} - \lambda_c^n}{\Delta t^n} = -\frac{1}{2} \left( |u^{n+1}| + \alpha \left| \frac{e_c^{n+1} - e_c^n}{\Delta t^n} \right| \right) \lambda_c^{n+\frac{1}{2}} + \sqrt{K_c^{n+1}} \frac{e_c^{n+1} - e_c^n}{\Delta t^n} $$

$$ + \frac{n_0(e_c^n)}{\sqrt{K_c^{n+1}}} \left( T_0 - \frac{K_0 \lambda_c^{n+\frac{1}{2}}}{2\sqrt{K_c^{n+1}}} \right) |u^{n+1}|_+. \quad (35) $$