VanOsta2024

Model summary

Cheatsheet

Tutorials

API reference

The VanOsta2024 model was configured and first published by Van Osta et al.[1], building upon the CircAdapt model published by Walmsley et al.[2]. This rapid biomechanical closed-loop model focused on simulating the intricacies of the heart with a realistic pre- and afterload simulated by a simple closed-loop circulation. Refer to Fig. 1 for a schematic overview of the model.

The foundation of this physiology-driven approach is the one-fiber model [3]. This model establishes relationships between cavity pressure and volume to myocardial fiber stress and strain by assuming homogeneity of fiber stress and strain throughout the cardiac walls. The most important assumptions are the thick-walled rotationally symmetric geometry for cardiac walls and conservation of energy. More details can be found here.

Fig. 1 Schematic overview of the VanOsta2024 model.

Heart structure and function

The VanOsta2024 model contains various integrated modules representing sarcomere mechanics, cardiac chambers, cardiac valves, pericardium, major blood vessels, and systemic and pulmonary resistance.

Sarcomere mechanics

Sarcomere mechanics are phenomenologically described in the VanOsta2024 model using a modified three-element Hill contraction model [4]. In this sarcomere model, myofiber stress is a summation of active stress, generated by the sarcomere, and passive stress resulting from stretching of elastic structures, primarily attributed to the extracellular matrix (ECM) and titin. The model reproduces basic properties of length-, and time-dependent active stress generation in cardiac tissue. The fiber stress is determined by the rise of contractility (representing density of actin-myosin cross bridge formation) and the fiber strain. More detailed information on sarcomere mechanics can be found in the Patch module.

Atrial representation

Each atrium is simulated by a Chamber module, a spherical Cavity surrounded by a spherical Wall representing a cardiac chamber having no mechanical interaction with other chambers. Although mechanical interaction of the atrial walls with each other and with the ventricular walls occurs in the human heart, these interactions are relatively insignificant [5] and therefore neglected in the model.

Ventricular representation

The ventricles are simulated using the TriSeg module [6], where a pair of cardiac chambers share the septal wall. The ventricles are more complex since the right and left ventricle interact through their shared wall, the intraventricular septum, and the TriSeg module takes this interaction into account. In the Triseg module, the two ventricles are encapsulated by three Walls. These walls are in mechanical equilibrium in the hinge points, resulting in strong mechanical interaction.

Multipatch

In many examples of cardiac pathology, mechanical properties are inhomogeneously distributed across the various walls. To simulate such inhomogeneities, the Wall can be subdivided into multiple patches, each with specific mechanical properties [2] and equal wall curvature. Each patch is modeled by the Patch module, allowing for specific timing of mechanical activation and myocardial fiber properties.

Cardiac valves

The Valve module is used to model any orifice in the circulation where blood passes through. The passing blood flow is caused by a pressure drop over the valve, while incorporating inertia and Bernoulli’s principle. This module is used for the atrio-ventricular valves, ventricular-arterial valves, veno-atrial inlets, systemic and intracardiac shunts.

Pericardium

The non-linear constraining effect of the pericardium on cardiac cavities is simulated using the Bag module, adding an external pressure to all cardiac chambers.

Circulatory system

A reduced-order closed-loop circulation is implemented in the VanOsta2024 model, consisting of the systemic and pulmonary circulation. Large arteries and veins are lumped in a single Tube0D, which describes the non-linear pressure-area relation of large blood vessels. Blood flow through the microcirculation is simulated using the ArtVen module, which describes flow, dependent on a pressure drop and a non-linear peripheral resistance.

Model control systems

Timing control

The Timings module controls the timing of contraction during a physiological heartbeat. This module manages the atrial-ventricular delay and heartrate-dependent parameters.

Homeostatic regulation

The Pressure Flow Control module represents homeostatic pressure and flow regulation. It regulates mean arterial pressure and venous return (cardiac output) by adjusting peripheral resistance and total blood volume.

Adaptation

The CircAdapt model allows for simulation of structural responses to chronic load by using adaptation rules [7][8]. This adaptation is optional and is recommended to use when creating a healthy reference parameterization or when simulating structural remodeling in compensated chronic pathophysiological situations.

Default parameterization

The model is initiated with a default parameterization which represents healthy physiology (Cheatsheet). This can be visualized using the default plot statement as shown below and can be generated using the following code:

import circadapt
import matplotlib.pyplot as plt

model = circadapt.VanOsta2024()
model.run(stable=True)
model.plot()
plt.show()

(Source code, png, hires.png, pdf)

Fig. 2 Default visualization of the model with its reference parameterization.

This default parameterization gives the following scalar output:

Table 1 Model compared to reference values

	Model	Literature	Unit	Source
LVEDV	120.5	121 pm 31	mL	Addetia et al.[9]
LVESV	48.3	47 pm 15	mL	Addetia et al.[9]
LVEF	59.9	61 pm 5	%	Addetia et al.[9]
RVEDV	109.2	91 [61 - 150]	mL	Maffessanti et al.[10]
RVESV	37.0	35 [16 - 72]	mL	Maffessanti et al.[10]
RVEF	66.1	62 [47 - 77]	%	Maffessanti et al.[10]
SBP	119.1	121 pm 12	mmHg	Addetia et al.[9]
DBP	70.3	74 pm 9	mmHg	Addetia et al.[9]
mPAP	17.2	[8 - 20]	mmHg	Humbert et al.[11]
mLAP	7.2	0 - 15	mmHg	Humbert et al.[11]
mRAP	2.0	2 - 6	mmHg	Humbert et al.[11]
AVdelay	150.0		ms

Model assumptions and parameterizations

Todo

• Present the initial parameterization.

• Discuss the model verification process.

• Explain adaptation protocol

Solving strategy

Each object in this model is derived from one of the CircAdapt modules. Together these modules define a closed loop system of ordinary differential equations that describe cardiovascular mechanics and hemodynamics. The system is solved using one of the implemented solvers. By default, an Adams Moulton solver is applied with a time step size of 0.001 ms.

The modular architecture is illustrated in Fig. 4. Four levels of organization are represented:

• Cavity modules describe pressures as a function of cavity volume and flow, including vascular elements such as Tube0D for wave transmission.

• Wall modules represent cardiac or vascular walls, where geometry and stress are linked through Laplace relations and material properties.

• Patch modules capture regional mechanics within a wall, with active and passive stresses depending on strain and activation.

• Cell modules describe excitation contraction coupling, relating calcium dynamics to active stress generation.

Information flows between modules: cavity pressure depends on wall tension, wall stress depends on patch stress, and patch stress is generated from cell level calcium dynamics. This hierarchical coupling ensures that organ level behavior emerges from interactions across scales, from cell to circulation.

Fig. 4 Overview of the hierarchical organization of CircAdapt modules.

References

[1]

Nick Van Osta, Gitte Van Den Acker, Tim Van Loon, Theo Arts, Tammo Delhaas, and Joost Lumens. Numerical accuracy of closed-loop steady state in a zero-dimensional cardiovascular model. Philosophical Transactions A, 383(2293):20240208, 2025.

[2] (1,2)

John Walmsley, Theo Arts, Nicolas Derval, Pierre Bordachar, Hubert Cochet, Sylvain Ploux, Frits W. Prinzen, Tammo Delhaas, and Joost Lumens. Fast simulation of mechanical heterogeneity in the electrically asynchronous heart using the multipatch module. PLOS Computational Biology, 11(7):1–23, 07 2015.

[3]

Theo Arts, PH Bovendeerd, Frits W Prinzen, and Robert 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, 1991.

[4]

Y.C. Fung. Biomechanics, Mechanical Properties of Living Tissues. Springer-Verlag, New York, second edition, 1993.

[5]

James A Goldstein, James S Tweddell, Benico Barzilai, Yoko Yagi, Allan S Jaffe, and James L Cox. Hemodynamic effects of atrial interaction. Coronary Artery Disease, 4(6):545–553, 1993.

[6]

Joost Lumens, Tammo Delhaas, Borut Kirn, and Theo Arts. Three-wall segment (triseg) model describing mechanics and hemodynamics of ventricular interaction. Annals of biomedical engineering, 37:2234–2255, 2009.

[7]

Theo Arts, Tammo Delhaas, Peter Bovendeerd, Xander Verbeek, and Frits W Prinzen. Adaptation to mechanical load determines shape and properties of heart and circulation: the circadapt model. American Journal of Physiology-Heart and Circulatory Physiology, 288(4):H1943–H1954, 2005.

[8]

Theo Arts, Joost Lumens, Wilco Kroon, and Tammo Delhaas. Control of whole heart geometry by intramyocardial mechano-feedback: a model study. PLoS computational biology, 8(2):e1002369, 2012.

[9] (1,2,3,4,5)

Karima Addetia, Tatsuya Miyoshi, Vivekanandan Amuthan, Rodolfo Citro, Masao Daimon, Pedro Gutierrez Fajardo, Ravi R Kasliwal, James N Kirkpatrick, Mark J Monaghan, Denisa Muraru, and others. Normal values of left ventricular size and function on three-dimensional echocardiography: results of the world alliance societies of echocardiography study. Journal of the American Society of Echocardiography, 35(5):449–459, 2022.

[10] (1,2,3)

Francesco Maffessanti, Denisa Muraru, Roberta Esposito, Paola Gripari, Davide Ermacora, Ciro Santoro, Gloria Tamborini, Maurizio Galderisi, Mauro Pepi, and Luigi P Badano. Age-, body size-, and sex-specific reference values for right ventricular volumes and ejection fraction by three-dimensional echocardiography: a multicenter echocardiographic study in 507 healthy volunteers. Circulation: Cardiovascular Imaging, 6(5):700–710, 2013.

[11] (1,2,3)

Marc Humbert, Gabor Kovacs, Marius M Hoeper, Roberto Badagliacca, Rolf MF Berger, Margarita Brida, Jørn Carlsen, Andrew JS Coats, Pilar Escribano-Subias, Pisana Ferrari, and others. 2022 esc/ers guidelines for the diagnosis and treatment of pulmonary hypertension: developed by the task force for the diagnosis and treatment of pulmonary hypertension of the european society of cardiology (esc) and the european respiratory society (ers). endorsed by the international society for heart and lung transplantation (ishlt) and the european reference network on rare respiratory diseases (ern-lung). European heart journal, 43(38):3618–3731, 2022.

Python implementation

class circadapt.VanOsta2024(solver: str = None, path_to_circadapt: str = None, model_state: dict = None)

Bases: Model, ModelAdapt

adapt(options: dict = {}, output: bool = False, verbose: bool = False) → None

Run adaptation protocol.

The adaptation protocol runs n_cycles cycles. Each cycle has two phases, namely rest adaptation and exercise adaptation. First, in exercise, the Patches and vessel wall volumes are adapted to load. In rest, vessels are adapted to flow.

Parameters

optionsdictionary, optional

Options for the protocol. To set options, use the function get_adapt_options(), which is used by default. The default is {}.

outputbool, optional

DESCRIPTION. The default is False.

Returns

senses_resultsTYPE

Only when output=True.

senses_normTYPE

Only when output=True.

actor_resultsTYPE

Only when output=True.

actor_vwallTYPE

Only when output=True.

adapt_exercise(options)

Trigger all excercise adaptation functions.

adapt_rest(options)

Trigger all rest adaptation functions.

add(comp_type: str) → bool

add_component(comp_type: str, comp_name: str, base: str = '') → bool

Add component to CircAdapt object.

Parameters

comp_type: str

Type of object to create in the ComponentFactory

comp_name: str

Name of the new object to create

base: char, optional (default=’’)

Parent object of new component

Returns

is_success (bool)

add_smart_component(smart_component, base=None, **kwargs)

assert_1_beat()

Test status after 1 beat, used in unit tests.

build()

Build model.

Add and link components in the model. By default, the model is empty. This function will be used in derived model builts.

build_artven(base=None, **kwargs)

build_heart(base=None, options=None)

build_pfc(base=None, options=None)

build_timings(base=None, options=None)

calculate_and_set_matrix(verbose=False)

check_build()

copy()

Return a copy of itself.

get(*arg, **kwarg)

get_adapt_options()

Get default adapt options.

get_matrix(verbose=False)

get_model_reference(model_state: any = None) → dict

Return reference model state as a dictionary.

If model_state is string, it is assumed to be a file location. This file will be loaded and returned. If is none, the default reference will be loaded. If no default reference available, it will be created and stored in the current active directory.

Parameters

model_statestr or None, optional

Location of model state. The default is None.

Raises

ValueError

DESCRIPTION.

Returns

dict

Dictonary with model state.

get_unittest_results(model)

Real-time results after initializing and running 1 beat.

get_unittest_targets()

Hardcoded results after initializing and running 1 beat.

health_check()

Manually perform tests to check compatibility of model with version.

install()

is_stable() → bool

Check if simulation is hemodynamically stable.

is_success() → bool

Return false if vectors contain nan.

Returns

bool

load(filename: str) → None

Load a model dataset from a filename.

Parameters

filename: str

Path to file that will be loaded. The extension must be .npy or .mat.

load_plugin_components(path_to_library: str)

Load a plugin

Parameters

path_to_librarystr

Path to library.

load_reference()

Load reference of current model from the package.

model_export(style: str = None) → dict

Return the stored model state.

Parameters

style: str (optional)

Style to follow. If not given, the default style from the model is used.

Returns

dict

model_import(model_state, check_model_state=False) → None

Load model_state into CircAdapt object.

Style and model_state version is automatically recognized.

Parameters

model_state: dict

model_state with data to set

obj: str

Object to set

plot(fig=None)

Simple plot to illustrate the model state.

plot_extended(fig=None)

Extended plot to illustrate the model state.

run(*arg, **kwarg)

save(filename: str, ext: str = None) → None

Save model to filename with extention.

Parameters

filename: str

Filename to save file to. Filename must have extention .npy or .mat

ext: str (optional)

Extention of the file. If not given, the last 3 letters of the filename is used to determine the save method.

set(*arg, **kwarg)

set_component(par, *arg, **kwarg)

set_reference()

stepper(t)

trigger(par) → bool

Trigger a function.

Objects in the model may have a function that can be triggered. This can be done by triggering the function using the ‘par’ parameter similar to the set function, only without a parameter.

Parameters

parstr

Function in object that will be triggered. It has a similar form to the set/get par, i.e. ‘Component.subcomponent.function_name’

Returns

is_success (bool)

Source code

class VanOsta2024(Model, ModelAdapt):
    def __init__(self,
                 solver: str = None,
                 path_to_circadapt: str = None,
                 model_state: dict = None,
                 ):
        if solver is None:
            solver = 'adams_moulton'

        self._local_save_reference = True

        ModelAdapt.__init__(self)
        Model.__init__(self,
                       solver,
                       path_to_circadapt=path_to_circadapt,
                       model_state=model_state,
                       )

    def _update_settings(self):
        # rename model components for easy input/output
        self._module_rename['PressureFlowControl'] = 'PFC'

    def build(self):
        pass
        # Circulation
        self.add_smart_component('ArtVen', build='SystemicCirculation')
        self.add_smart_component('ArtVen', build='PulmonaryCirculation')
        self.add_smart_component('Heart', patch_type='Patch', valve_type='Valve')
        self.add_smart_component('Timings')
        self.add_smart_component('PressureFlowControl')

        self.set_component('PFC.circulation_volume_object', 'SyArt')
        self.set_component('PFC.circulation_volume_object', 'SyVen')
        self.set_component('PFC.circulation_volume_object', 'PuArt')
        self.set_component('PFC.circulation_volume_object', 'PuVen')
        self.set_component('PFC.circulation_volume_object', 'Peri.La')
        self.set_component('PFC.circulation_volume_object', 'Peri.Ra')
        self.set_component('PFC.circulation_volume_object', 'Peri.TriSeg.cLv')
        self.set_component('PFC.circulation_volume_object', 'Peri.TriSeg.cRv')

        # # manually set papillary muscles
        self.set_component("Peri.RaRv.wPapMus", "Peri.TriSeg.wRv")
        self.set_component("Peri.LaLv.wPapMus", "Peri.TriSeg.wLv")

    def set_reference(self):
        self['Chamber']['buckling'] = True
        self['Valve']['soft_closure'] = True
        self['Valve']['papillary_muscles'] = False

        self['Solver']['dt'] = 0.001
        self['Solver']['dt_export'] = 0.002
        self.set('Solver.order',  2)

        self.set('Model.t_cycle', 0.85)
        self['ArtVen']['p0'] = np.array([6306.25832487, 1000.        ])
        self['ArtVen']['q0'] = np.array([4.5e-05, 4.5e-05])
        self['ArtVen']['k'] = np.array([1., 2.])
        self['Tube0D']['l'] = np.array([0.4, 0.4, 0.2, 0.2])
        self['Tube0D']['A_wall'] = np.array([1.12362733e-04, 6.57883944e-05, 9.45910889e-05, 8.22655361e-05])
        self['Tube0D']['k'] = np.array([1.66666667, 2.33333333, 1.66666667, 2.33333333])
        self['Tube0D']['p0'] = np.array([12162.50457811,   287.7083132 ,  2132.51755623,   830.54673184])
        self['Tube0D']['A0'] = np.array([0.0004983 , 0.00049909, 0.00047138, 0.00050803])
        self['Tube0D']['target_wall_stress'] = np.array([500000., 500000., 500000., 500000.])
        self['Tube0D']['target_mean_flow'] = np.array([0.17, 0.17, 0.17, 0.17])
        self['Bag']['k'] = np.array([10.])
        self['Bag']['p_ref'] = np.array([1000.])
        self['Bag']['V_ref'] = np.array([0.00054267])
        self['Patch']['Am_ref'] = np.array([0.00425687, 0.00401573, 0.00966859, 0.00289936, 0.01084227])
        self['Patch']['V_wall'] = np.array([4.46069398e-06, 2.14548521e-06, 7.35720515e-05, 1.88904978e-05, 3.67720116e-05,])
        self['Patch']['v_max'] = np.array([14., 14.,  7.,  7.,  7.])
        self['Patch']['l_se0'] = np.array([0.04, 0.04, 0.04, 0.04, 0.04])
        self['Patch']['l_s0'] = np.array([1.8, 1.8, 1.8, 1.8, 1.8])
        self['Patch']['l_s_ref'] = np.array([2., 2., 2., 2., 2.])
        self['Patch']['dl_s_pas'] = np.array([0.6, 0.6, 0.6, 0.6, 0.6])
        self['Patch']['Sf_pas'] = np.array([2248.53598   , 2684.76100348,  731.24545453,  729.06063771, 749.47694522,])
        self['Patch']['tr'] = np.array([0.4 , 0.4 , 0.25, 0.25, 0.25])
        self['Patch']['td'] = np.array([0.4 , 0.4 , 0.25, 0.25, 0.25])
        self['Patch']['time_act'] = np.array([0.15 , 0.15 , 0.425, 0.425, 0.425])
        self['Patch']['Sf_act'] = np.array([ 84000.,  84000., 120000., 120000., 120000.])
        self['Patch']['fac_Sf_tit'] = np.array([0.01, 0.01, 0.01, 0.01, 0.01])
        self['Patch']['k1'] = np.array([10., 10., 10., 10., 10.])
        self['Patch']['dt'] = np.array([0., 0., 0., 0., 0.])
        self['Patch']['C_rest'] = np.array([0., 0., 0., 0., 0.])
        self['Patch']['l_si0'] = np.array([1.51, 1.51, 1.51, 1.51, 1.51])
        self['Patch']['LDAD'] = np.array([1.057, 1.057, 1.057, 1.057, 1.057])
        self['Patch']['ADO'] = np.array([0.65, 0.65, 0.65, 0.65, 0.65])
        self['Patch']['LDCC'] = np.array([4., 4., 4., 4., 4.])
        self['Patch']['SfPasMaxT'] = np.array([320000., 320000.,  16000.,  16000.,  16000.])
        self['Patch']['SfPasActT'] = np.array([40000., 40000., 20000., 20000., 20000.])
        self['Patch']['FacSfActT'] = np.array([0.8, 0.8, 1. , 1. , 1. ])
        self['Patch']['LsPasActT'] = np.array([2.42, 2.42, 2.42, 2.42, 2.42])
        self['Patch']['adapt_gamma'] = np.array([0.5, 0.5, 0.5, 0.5, 0.5])
        self['Patch']['transmat00'] = np.array([-0.5751, -0.5751, -0.5751, -0.5751, -0.5751])
        self['Patch']['transmat01'] = np.array([-0.7851, -0.7851, -0.7851, -0.7851, -0.7851])
        self['Patch']['transmat02'] = np.array([0.6063, 0.6063, 0.6063, 0.6063, 0.6063])
        self['Patch']['transmat03'] = np.array([-0.5565, -0.5565, -0.5565, -0.5565, -0.5565])
        self['Patch']['transmat10'] = np.array([-0.1279, -0.1279, -0.1279, -0.1279, -0.1279])
        self['Patch']['transmat11'] = np.array([0.0999, 0.0999, 0.0999, 0.0999, 0.0999])
        self['Patch']['transmat12'] = np.array([0.2066, 0.2066, 0.2066, 0.2066, 0.2066])
        self['Patch']['transmat13'] = np.array([-1.8441, -1.8441, -1.8441, -1.8441, -1.8441])
        self['Patch']['transmat20'] = np.array([-0.1865, -0.1865, -0.1865, -0.1865, -0.1865])
        self['Patch']['transmat21'] = np.array([-0.201, -0.201, -0.201, -0.201, -0.201])
        self['Patch']['transmat22'] = np.array([1.3195, 1.3195, 1.3195, 1.3195, 1.3195])
        self['Patch']['transmat23'] = np.array([-11.8745, -11.8745, -11.8745, -11.8745, -11.8745])
        self['Valve']['adaptation_A_open_fac'] = np.array([1., 1., 1., 1., 1., 1.])
        self['Valve']['A_open'] = np.array([0.00049916, 0.00047155, 0.00047155, 0.00050805, 0.00049835,
               0.00049835])
        self['Valve']['A_leak'] = np.array([0.00049916, 1.e-09, 1.e-09, 0.00050805, 1.e-09, 1.e-09])
        self['Valve']['l'] = np.array([0.01260512, 0.01225144, 0.01225144, 0.01271679, 0.01259479,
               0.01259479])
        self['Valve']['L_fac_prox'] = np.array([0.75, 0.75, 0.75, 0.75, 0.75, 0.75])
        self['Valve']['L_fac_dist'] = np.array([0.75, 0.75, 0.75, 0.75, 0.75, 0.75])
        self['Valve']['L_fac_valve'] = np.array([1.5, 1.5, 1.5, 1.5, 1.5, 1.5])
        self['Valve']['rho_b'] = np.array([1050., 1050., 1050., 1050., 1050., 1050.])
        self['Valve']['papillary_muscles'] = np.array([False, False, False, False, False, False])
        self['Valve']['papillary_muscles_slope'] = np.array([100., 100., 100., 100., 100., 100.])
        self['Valve']['papillary_muscles_min'] = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
        self['Valve']['papillary_muscles_A_open_fac'] = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
        self['Valve']['soft_closure'] = np.array([ True,  True,  True,  True,  True,  True])
        self['Valve']['fraction_A_open_Aext'] = np.array([0.9, 0.9, 0.9, 0.9, 0.9, 0.9])
        self['Timings']['time_fac'] = np.array([1.])
        self['Timings']['tau_av'] = np.array([0.150025])
        self['Timings']['dtau_av'] = np.array([0.])
        self['Timings']['law_tau_av'] = np.array([1])
        self['Timings']['law_Ra2La'] = np.array([1])
        self['Timings']['law_ta'] = np.array([1])
        self['Timings']['law_tv'] = np.array([1])
        self['Timings']['c_tau_av0'] = np.array([0.])
        self['Timings']['c_tau_av1'] = np.array([0.1765])
        self['Timings']['c_ta_rest'] = np.array([0.])
        self['Timings']['c_ta_tcycle'] = np.array([0.17647059])
        self['Timings']['c_tv_rest'] = np.array([0.1])
        self['Timings']['c_tv_tcycle'] = np.array([0.4])
        self['PFC']['p0'] = np.array([12200.])
        self['PFC']['q0'] = np.array([8.5e-05])
        self['PFC']['stable_threshold'] = np.array([0.0001])
        self['PFC']['is_active'] = np.array([ True])
        self['PFC']['fac'] = np.array([1.])
        self['PFC']['alpha'] = np.array([1.])
        self['PFC']['epsilon'] = np.array([0.4])
        self['TriSeg']['tau'] = np.array([2.])
        self['TriSeg']['ratio_septal_LV_Am'] = np.array([0.3])
        self['TriSeg']['max_number_of_iterations'] = np.array([100.])
        self['TriSeg']['thresh_F'] = np.array([0.001])
        self['TriSeg']['thresh_dV'] = np.array([1.e-09])
        self['TriSeg']['thresh_dY'] = np.array([1.e-06])
        self['TriSeg']['collapse_eps'] = 0.0
        self['TriSeg']['collapse_start'] = 0.25

        self.run(stable=True)
        # self.adapt()

        # set target volume
        self['PFC']['target_volume'] = self['PFC']['circulation_volume']

        return


    def get_unittest_targets(self):
        """Hardcoded results after initializing and running 1 beat."""
        return {
            'LVEDV': 120.3,
            'LVESV':  48.0,
            }

    def get_unittest_results(self, model):
        """Real-time results after initializing and running 1 beat."""
        LVEDV = np.max(model['Cavity']['V'][:, 'cLv'])*1e6
        LVESV = np.min(model['Cavity']['V'][:, 'cLv'])*1e6
        return {
            'LVEDV': LVEDV,
            'LVESV': LVESV,
            }


    def plot(self, fig=None):
        # TODO
        self.plot_extended(fig)

    def plot_extended(self, fig=None):
        if len(self['Solver']['t'])==0:
            raise RuntimeError('Can not plot the model, as no data is available. Did you simulate the model?')

        if fig is None:
            fig = 1
        if isinstance(fig, int):
            fig = plt.figure(fig, clear=True, figsize=(12, 8))

        # Settings
        grid_size = [32, 32]

        def get_lim(module, signal, locs=slice(None, None, None)):
            signal = self[module][signal][:, locs]
            lim = np.array([np.min(signal), np.max(signal)])
            lim += np.array([-1, 1]) * 0.1*np.diff(lim)
            return lim

        lim_V = get_lim('Cavity', 'V', ['cRv', 'Ra', 'La', 'cLv']) * 1e6
        lim_V[0] = 0
        lim_p = get_lim('Cavity', 'p', ['cRv', 'Ra', 'La', 'cLv']) / 133
        lim_p[0] = np.min([lim_p[0], 0])
        lim_Ls = get_lim('Patch', 'l_s')
        lim_Sf = get_lim('Patch', 'Sf') * 1e-3
        lim_q = get_lim('Valve', 'q',
                        ['LaLv', 'RaRv', 'LvSyArt', 'RaPuArt']) * 1e6

        all_lim = [lim_V, lim_p, lim_Ls, lim_Sf, lim_q]
        if (np.any(np.isnan(all_lim)) or np.any(np.isinf(all_lim))):
            lim_V = [0, 200]
            lim_p = [0, 150]
            lim_Ls = [1.5, 2.0]
            lim_Sf = [0, 100]
            lim_q = [-1e-3, 1e-3]

        # Pressure Volume plot
        axPV = plt.subplot2grid(grid_size, (0, 17), rowspan=15, colspan=15, fig=fig)
        axPV.plot(self['Cavity']['V'][:, 'cLv']*1e6,
                  self['Cavity']['p'][:, 'cLv']/133,
                  c=colors['Lv'],
                  )
        axPV.plot(self['Cavity']['V'][:, 'cRv']*1e6,
                  self['Cavity']['p'][:, 'cRv']/133,
                  c=colors['Rv'],
                  )
        axPV.plot(self['Cavity']['V'][:, 'La']*1e6,
                  self['Cavity']['p'][:, 'La']/133,
                  c=colors['La'],
                  )
        axPV.plot(self['Cavity']['V'][:, 'Ra']*1e6,
                  self['Cavity']['p'][:, 'Ra']/133,
                  c=colors['Ra'],
                  )
        axPV.spines[['top', 'right']].set_visible(False)
        axPV.set_title('Pressure-Volume loop', weight='bold')
        axPV.set_xlabel('Volume [mL]')
        axPV.set_ylabel('Pressure [mmHg]')
        axPV.spines[['bottom', 'left']].set_position(('outward', 5))

        ylabel_x_left = -0.25
        ylabel_x_right = 1.25

        # Volumes
        t = self['Solver']['t']*1e3

        axVRv = plt.subplot2grid(grid_size, (0, 0), rowspan=8, colspan=6, fig=fig)
        axVRv.plot(t, self['Cavity']['V'][:, 'cRv']*1e6,
                   c=colors['Rv'],
                   )
        axVRv.plot(t, self['Cavity']['V'][:, 'Ra']*1e6,
                   c=colors['Ra'],
                   )
        axVRv.set_ylim(lim_V)
        axVRv.set_ylabel('Volume\n[mL]')
        axVRv.spines[['top', 'right']].set_visible(False)
        axVRv.set_title('Right Heart', weight='bold')
        # axVRv.set_xticks([])
        axVRv.tick_params(axis='both', direction='in')
        axVRv.yaxis.set_label_coords(ylabel_x_left, 0.5)


        axVLv = plt.subplot2grid(grid_size, (0, 6), rowspan=8, colspan=6, fig=fig)
        axVLv.plot(t, self['Cavity']['V'][:, 'cLv']*1e6,
                   c=colors['Lv'],
                   )
        axVLv.plot(t, self['Cavity']['V'][:, 'La']*1e6,
                   c=colors['La'],
                   )
        axVLv.set_ylabel('Volume\n[mL]')
        axVLv.set_ylim(lim_V)
        axVLv.yaxis.set_ticks_position('right')
        axVLv.yaxis.set_label_position('right')
        axVLv.spines['right'].set_position(('outward', 0))
        axVLv.spines[['top', 'left']].set_visible(False)
        axVLv.set_title('Left Heart', weight='bold')
        # axVLv.set_xticks([])
        axVLv.tick_params(axis='both', direction='in')
        axVLv.yaxis.set_label_coords(ylabel_x_right, 0.5)

        # Pressures
        axpRv = plt.subplot2grid(
            grid_size, (8, 0), rowspan=8, colspan=6, fig=fig, sharex=axVRv)
        axpRv.plot(t, self['Cavity']['p'][:, 'cRv']/133,
                   c=colors['Rv'],
                   )
        axpRv.plot(t, self['Cavity']['p'][:, 'Ra']/133,
                   c=colors['Ra'],
                   )
        axpRv.plot(t, self['Cavity']['p'][:, 'PuArt']/133,
                   c=colors['PuArt'],
                   )
        axpRv.spines[['top', 'right']].set_visible(False)
        # axpRv.set_xticks([])
        axpRv.tick_params(axis='both', direction='in')
        axpRv.set_ylim(lim_p)
        axpRv.set_ylabel('Pressure\n[mmHg]')
        axpRv.yaxis.set_label_coords(ylabel_x_left, 0.5)

        axpLv = plt.subplot2grid(grid_size, (8, 6), rowspan=8, colspan=6, fig=fig, sharex=axVLv)
        axpLv.plot(t, self['Cavity']['p'][:, 'cLv']/133,
                   c=colors['Lv'],
                   )
        axpLv.plot(t, self['Cavity']['p'][:, 'La']/133,
                   c=colors['La'],
                   )
        axpLv.plot(t, self['Cavity']['p'][:, 'SyArt']/133,
                   c=colors['SyArt'],
                   )
        axpLv.yaxis.set_ticks_position('right')
        axpLv.yaxis.set_label_position('right')
        axpLv.spines['right'].set_position(('outward', 0))
        axpLv.spines[['top', 'left']].set_visible(False)
        # axpLv.set_xticks([])
        axpLv.tick_params(axis='both', direction='in')
        axpLv.set_ylim(lim_p)
        axpLv.set_ylabel('Pressure\n[mmHg]')
        axpLv.yaxis.set_label_coords(ylabel_x_right, 0.5)

        # Valves
        ax = plt.subplot2grid(grid_size, (16, 0), rowspan=6, colspan=6, fig=fig, sharex=axVRv)
        ax.plot(t, self['Valve']['q'][:, 'RaRv']*1e6,
                c=colors['RaRv'],
                )
        ax.plot(t, self['Valve']['q'][:, 'RvPuArt']*1e6,
                c=colors['RvPuArt'],
                )
        ax.spines[['top', 'right']].set_visible(False)
        ax.set_ylim(lim_q)
        # ax.set_xticks([])
        ax.set_ylabel('Flow\n[mL/s]')
        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)

        ax = plt.subplot2grid(grid_size, (16, 6), rowspan=6, colspan=6, fig=fig, sharex=axVLv)
        ax.plot(t, self['Valve']['q'][:, 'LaLv']*1e6,
                   c=colors['LaLv'],
                   )
        ax.plot(t, self['Valve']['q'][:, 'LvSyArt']*1e6,
                   c=colors['LvSyArt'],
                   )
        ax.spines[['top', 'left']].set_visible(False)
        ax.set_ylim(lim_q)
        ax.yaxis.set_ticks_position('right')
        ax.yaxis.set_label_position('right')
        # ax.set_xticks([])
        ax.set_ylabel('Flow\n[mL/s]')
        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)

        # Stress
        ax = plt.subplot2grid(grid_size, (22, 0), rowspan=4, colspan=6, fig=fig, sharex=axVRv)
        ax.plot(t, self['Patch']['Sf'][:, 'pRv0']*1e-3,
                   c=colors['Rv'],
                   )
        ax.plot(t, self['Patch']['Sf'][:, 'pRa0']*1e-3,
                   c=colors['Ra'],
                   )
        ax.spines[['top', 'right']].set_visible(False)
        # ax.set_xticks([])
        ax.set_ylim(lim_Sf)
        ax.set_ylabel('Total\nstress [kPa]')
        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)

        ax = plt.subplot2grid(grid_size, (22, 6), rowspan=4, colspan=6, fig=fig, sharex=axVLv)
        ax.plot(t, self['Patch']['Sf'][:, 'pLv0']*1e-3,
                   c=colors['Lv'],
                   )
        ax.plot(t, self['Patch']['Sf'][:, 'pSv0']*1e-3,
                   c=colors['Sv'],
                   )
        ax.plot(t, self['Patch']['Sf'][:, 'pLa0']*1e-3,
                   c=colors['La'],
                   )
        ax.spines[['top', 'left']].set_visible(False)
        ax.yaxis.set_ticks_position('right')
        ax.yaxis.set_label_position('right')
        # ax.set_xticks([])
        ax.set_ylim(lim_Sf)
        ax.set_ylabel('Total\nstress [kPa]')
        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)

        # Sarcomere Length
        ax = plt.subplot2grid(grid_size, (26, 0), rowspan=6, colspan=6, fig=fig, sharex=axVRv)
        ax.plot(t, self['Patch']['l_s'][:, 'pRv0'],
                   c=colors['Rv'],
                   )
        ax.plot(t, self['Patch']['l_s'][:, 'pRa0'],
                   c=colors['Ra'],
                   )
        ax.spines[['top', 'right']].set_visible(False)
        ax.spines['bottom'].set_position(('outward', 5))
        ax.set_ylim(lim_Ls)
        ax.set_ylabel('Sarcomere\nlength [$\\mu$m]')
        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)

        ax = plt.subplot2grid(grid_size, (26, 6), rowspan=6, colspan=6, fig=fig, sharex=axVLv)
        ax.plot(t, self['Patch']['l_s'][:, 'pLv0'],
                   c=colors['Lv'],
                   )
        ax.plot(t, self['Patch']['l_s'][:, 'pSv0'],
                   c=colors['Sv'],
                   )
        ax.plot(t, self['Patch']['l_s'][:, 'pLa0'],
                   c=colors['La'],
                   )
        ax.spines[['top', 'left']].set_visible(False)
        ax.spines['bottom'].set_position(('outward', 5))
        ax.yaxis.set_ticks_position('right')
        ax.yaxis.set_label_position('right')
        ax.set_ylim(lim_Ls)
        ax.set_ylabel('Sarcomere\nlength [$\\mu$m]')
        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)

        # ax.set_xlabel('Time [ms]')
        # ax.xaxis.set_label_coords(0, -0.3)


        # Plot TriSeg
        titles = ['Pre-A', 'Onset QRS', 'Peak LV \n pressure', 'AV close']
        idx = [0,
               np.argmax(np.diff(self['Patch']['C'][:, 'pLv0'])>0),
               np.argmax(self['Cavity']['p'][:, 'cLv']),
               len(t) - 1 - np.argmax(
                   np.diff(self['Valve']['q'][:, 'LvSyArt'][::-1])>0)
               ]

        for i in range(4):
            ax = plt.subplot2grid(grid_size, (26, 16+4*i),
                                  rowspan=5, colspan=4, fig=fig)
            triseg2022(self, ax, idx[i],
                       colors=[colors['Lv'], colors['Sv'], colors['Rv']])
            ax.spines[['top', 'right', 'bottom', 'left']].set_visible(False)
            ax.set_xticks([])
            ax.set_yticks([])
            plt.xlabel(titles[i], fontsize=12)

        # Plot settings
        plt.subplots_adjust(
            top=0.96,
            bottom=0.05,
            left=0.075,
            right=0.98,
            hspace=5,
            wspace=0.5)
        plt.draw()

        # Stress strain
        ax_right_stress_strain = plt.subplot2grid(grid_size, (18, 17),
                              rowspan=7, colspan=7, fig=fig)
        ax_left_stress_strain = plt.subplot2grid(grid_size, (18, 24),
                              rowspan=7, colspan=7, fig=fig)


        ax_right_stress_strain.plot(
            self['Patch']['Ef'][:, 'pRa0'],
            self['Patch']['Sf'][:, 'pRa0']*1e-3,
            c=colors['Ra'],
            )
        ax_right_stress_strain.plot(
            self['Patch']['Ef'][:, 'pRv0'],
            self['Patch']['Sf'][:, 'pRv0']*1e-3,
            c=colors['Rv'],
            )
        ax_left_stress_strain.plot(
            self['Patch']['Ef'][:, 'pLa0'],
            self['Patch']['Sf'][:, 'pLa0']*1e-3,
            c=colors['La'],
            )
        ax_left_stress_strain.plot(
            self['Patch']['Ef'][:, 'pLv0'],
            self['Patch']['Sf'][:, 'pLv0']*1e-3,
            c=colors['Lv'],
            )
        ax_left_stress_strain.plot(
            self['Patch']['Ef'][:, 'pSv0'],
            self['Patch']['Sf'][:, 'pSv0']*1e-3,
            c=colors['Sv'],
            )

        ylim = [np.min(self['Patch']['Sf'])*1e-3,
                np.max(self['Patch']['Sf'])*1e-3]
        ylim_ptp = np.ptp(ylim)*0.025
        ylim[0] -= ylim_ptp
        ylim[1] += ylim_ptp
        ax_right_stress_strain.set_ylim(ylim)
        ax_left_stress_strain.set_ylim(ylim)
        xlim = [np.min(self['Patch']['Ef']),
                np.max(self['Patch']['Ef'])]
        xlim_ptp = np.ptp(xlim)*0.025
        xlim[0] -= xlim_ptp
        xlim[1] += xlim_ptp
        ax_right_stress_strain.set_xlim(xlim)
        ax_left_stress_strain.set_xlim(xlim)

        ax_right_stress_strain.spines[['left', 'bottom']].set_position(('outward', 5))
        ax_left_stress_strain.spines[['right', 'bottom']].set_position(('outward', 5))

        ax_right_stress_strain.spines[['top', 'right']].set_visible(False)
        ax_right_stress_strain.set_ylabel('Stress\n[kPa]', rotation=0, va='bottom', ha='right')
        ax_right_stress_strain.yaxis.set_label_coords(0., 1.05)
        ax_right_stress_strain.set_xlabel('Strain [-]')
        ax_left_stress_strain.spines[['top', 'left']].set_visible(False)
        ax_left_stress_strain.yaxis.set_ticks_position('right')
        ax_left_stress_strain.yaxis.set_label_position('right')
        ax_left_stress_strain.set_ylabel('Stress\n[kPa]', rotation=0, va='bottom', ha='left')
        ax_left_stress_strain.yaxis.set_label_coords(1., 1.05)
        ax_left_stress_strain.set_xlabel('Strain [-]')