1 Introduction
Computational modelling of cardiac EM is increasingly proposed and pursued as an effective tool for investigating cardiac physiology and numerous pathological conditions Augustin et al. (2016); Niederer et al. (2019). Computer models have been used for quantitative analysis, since they have the potential to assess valuable diagnostic information and represent a strong predictive tool for analysing the patientspecific response to a given treatment DiazZuccarini et al. (2013). However, this is a very open and challenging field of research, since cardiac function relies on complex biophysical aspects that are strongly interlinked Chabiniok et al. (2016). Indeed, the muscle contraction is stimulated by electrical activation and strongly interacts with intraventricular blood and the circulatory system to transport nutrients and clear waste products. As a consequence, efficient computational tools are required that account for the complex physiology and multiphysics of the heart and can be practically implemented and integrated into the clinic CorralAcero et al. (2020). In this context, numerous mathematical models have been reported in the literature that aim to model the effect of the circulatory system on cardiac function. These range from simple lumped zerodimensional (0D) Windkesseltype models Westerhof et al. (2009); Segers et al. (2003, 2008), to more refined 1D models Formaggia et al. (2003); Reymond et al. (2011); Manganotti et al. (2021). While 0D models are suitable as global models of afterload, providing a simple way to evaluate global cardiovascular dynamics, they do not consider any effects due to pulse wave transmission and thus are unable to predict important markers such as pulse wave velocity. In contrast, 1D blood flow models enable an efficient representation of the effects of pulse wave propagation and reflection in the circulation. Consequently, 1D models may be preferred over 0D models when local vascular changes or distributed properties, e.g., vessel branching, tapering, stenoses, and their effect on central pressure waveforms are studied and when pulse wave transmission effects are under investigation Van de Vosse and Stergiopulos (2011). To date, most 3D computational models of the heart proposed in the literature consider lumped 0D models as boundary conditions to model the circulation SainteMarie et al. (2006); Kerckhoffs et al. (2007); Nordsletten et al. (2011); Augustin et al. (2016); Hirschvogel et al. (2017); Regazzoni et al. (2020) for ease of implementation and calibration. To the best of our knowledge, no studies propose an effective and stable method to couple a 3D cardiac model to a 1D description of the circulation, which can provide new insights into the physiology and analyse the effects of disrupted wave propagation on cardiovascular dynamics. This is clinically relevant, since the impact of increased stiffness and pulsatile load on the circulation and cardiac performance have been clinically documented not only for cardiovascular events, but also for left and right ventricle dysfunctions Stevens et al. (2012). In order to address this need, in this methodological paper we present a novel in silico model based on a multidimensional approach, namely a 3D EM model of cardiac function together with a 1D model of the human circulation. The building blocks of the coupled model and the main aspects of model parameterisation are thoroughly described in the text. The coupling strategy is inspired by recent methods proposed for 3D0D models SainteMarie et al. (2006); Gurev et al. (2015) and is based on the solution of a saddlepoint problem for the volume and pressure in the cardiac cavity. A personalised 3D EM model of the LV Augustin et al. (2016) coupled to a 1D circulatory model Alastruey et al. (2012) is built and employed in numerical test cases to demonstrate robustness and stability of the numerical scheme. The physiological response of the coupled model is investigated in different conditions affecting pulse wave transmission, such as aortic stiffening, aortic narrowing and complex vessel networks involving numerous bifurcations. We show the ability of the coupled 3D1D model to correctly replicate known physiological behaviours related to pulse wave propagation. The computational cost of our coupled 3D1D model is comparable to that of standard 3D0D models, suggesting the use of our 3D1D model in clinical applications where wave transmission effects are investigated.
2 Methodology
2.1 3D Electromechanical cardiac model
Cardiac tissue is modelled as a hyperelastic, anisotropic, nearly incompressible material endowed with a nonlinear stressstrain relationship Guccione et al. (1995). Stresses associated with active contraction are characterised as orthotropic, and it is assumed that full contractile force acts along the myocyte fibre orientation, whereas 40% contractile force is directed along the sheet orientation Genet et al. (2014). Active stress generation is represented based on a simplified phenomenological contractile model Niederer et al. (2011). Active stresses in this model are only lengthdependent, and dependence on fibre velocity is neglected. A recent reactioneikonal (RE) model Neic et al. (2017) is considered for the generation of the electrical activation sequences that trigger active stress generation in the tissue. The hybrid RE model is a combination of a standard reactiondiffusion model (based on the monodomain equation) and an eikonal model. Spatiotemporal discretisation of all PDEs and the solvers for the resulting systems of equations are based on the Cardiac Arrhythmia Research Package (CARPentry) Augustin et al. (2016); Neic et al. (2017); Vigmond et al. (2008), built upon extensions of the freely available openCARP EP framework (http://www.opencarp.org). Further information on the EM model of cardiac function can be retrieved in A. For numerical details on the finite element discretisation, we refer the reader to E.1.
2.2 Circulatory system
The circulatory system imposes a load on the heart and, therefore, affects its mechanical activity. However, the interaction between the heart and the circulatory system is bidirectional, i.e., the outflow of blood from a heart cavity and pressure in the cavity depend on the current state of the circulatory system, and pressure and flow in the cardiovascular system are determined by the current state of the cavity itself. The full physics of this interaction is most accurately posed as a fluidstructure interaction (FSI) problem where pressure, , and flow velocity, , of the fluid are the coupling variables Karabelas et al. (2018); Fernández et al. (2007); Gerbeau (2012). Any perturbation in blood flow velocity and pressure changes the state of deformation of the heart and attached vessels. As a result, this change in strain implies a change in stress within the myocardial muscle. Conversely, any change in strain or compliance of the attached vessel or the heart changes the pressure and flow of the fluid. Even though such a distributed PDEbased approach may accurately describe this interaction, at the level of the whole circulatory system it is hardly feasible, for computational and structural reasons. In fact, a multibeat simulation using a 3D FSI model solution would be inconvenient for clinical application. Moreover, the use of 3D models of haemodynamics would require the identification of a much larger number of parameters, that is not easily attainable within clinical constraints. A reducedorder approximation of the circulatory system is used instead in this work, which based on a 1D model for the circulatory system and is particularly suitable for simulating blood flow in the aorta and the major arteries. The parameters in the circulatory model are identified and constrained by imagingbased measurements of the heart and blood flow and invasive blood pressure measurements. In this framework the nonlinear PDE EM model of the heart is coupled to the 1D model of the circulatory system using the lumped hydrostatic pressure in the cardiac cavity and the flow out of the cavity into the circulatory system as coupling variables.
2.2.1 1D mathematical model of the blood flow circulation in the arterial system
The 1D equations of blood flow in the human circulation can be derived from the NavierStokes equations after assuming axisymmetric flow in a cylindrical tube with thin wall Formaggia et al. (2003). In this work, we only take into account the arterial circulation and we consider a Voigttype viscoelastic constitutive law for the arterial wall Alastruey et al. (2012). A threeelement Windkessel lumped parameter model is employed as a terminal boundary condition to simulate the effect of peripheral vessels on pulse wave propagation in larger arteries. The numerical discretisation of the resulting system of equations relies upon a FEM discretisation (discontinuous Galerkin) in space and finite difference (explicit secondorder AdamsBashforth) in time Alastruey et al. (2012). The numerical solution of the resulting hyperbolic system of equations is performed with the solver Nektar1D (http://haemod.uk/nektar). We refer the reader to B for more information on 1D blood flow modelling. The arterial network is connected to the left ventricle (LV) cavity through a model of aortic valve (AV) dynamics, based on Bernoulli’s equation Mynard et al. (2012). See C for further details on the valve dynamics model.
3 3D1D coupling
Following the approach for 3D0D coupling proposed in Augustin et al. (2021) and inspired by Gurev et al. (2011); Kerckhoffs et al. (2007), the coupling condition between an EM model of the heart cavities and a reducedorder circulatory model imposes that the volume change in each cavity (left ventricle LV, right ventricle RV, left atrium LA and right atrium RA) is balanced with the volume change in the attached circulatory system. For the sake of generality, we cast in what follows the general framework to couple a fourchamber heart model with a closedloop circulatory system. Let us consider that approximately, at a given timepoint, we have constant pressures in each cavity , with . Then, the following equations are the starting point of the method:
(1) 
where is the cavity volume computed with the deformation using Eq. (2), and is the volume as predicted by the blood flow model for the intracavitary pressure . We write
for the vector of up to
pressure unknowns. Note, however, that in a purely EM simulation framework the fluid domain is not explicitly modelled. Therefore, the cavitary volume is not discretised. Instead, only the surface enclosing the cavitary volume is known. Nonetheless, if we assume that the entire surface of the cavitary volume is available, including the faces representing the valves, then we can compute for each instant of the cardiac cycle the enclosed volume from the surface by means of the divergence theorem(2) 
Then, the approach used to evolve the full system of equations of nonlinear elasticity Eq.(1) together with the coupling condition Eq.(2) is based on the resolution of a saddlepoint problem in the variables , corresponding to the displacement and the pressure in the cavity:
(3) 
which is valid for all vector fields smooth enough and satisfying the given boundary conditions, test functions that are for the cavity and otherwise, the duality pairing , and cavities . Note that can be physically interpreted as the rate of internal mechanical work, whereas takes into account the contribution of pressure loads. Using a Galerkin discretisation and the Newton method, the problem to be solved at each Newton–Raphson step becomes a block system in order to find and such that:
(4) 
with the updates
(5) 
and the solution vectors and at the th Newton step. In particular, we can retrieve the expression of and from the equations of the circulatory model, whereas the term only depends on the contribution of the EM model of the heart. In more detail, is approximated by a discrete derivative (finite difference) of the volume with respect to the cavity pressure. We refer the reader to D and E for further information on the coupling strategy and finite element formulation, respectively.
4 Simulation setup
4.1 Numerical framework
The coupling scheme and a C++ circulatory system module have been embedded in CARPentry, based on the software Nektar1D, in order to preserve computational efficiency and strong scalability. A Schur complement approach (see F) is considered to cast the coupling problem in a pure displacement formulation, in order to use solver methods already implemented Augustin et al. (2016). A generalised minimal residual method (GMRES) with an relative error reduction of is employed. The library PETSc Balay et al. (2018) and the incorporated solver suite hypre/BoomerAMG Henson and Yang (2002) are used for efficient preconditioning. A generalised scheme is considered for time integration, with spectral radius and damping parameters , , see Augustin et al. (2021) for further information.
4.2 Model parameterisation
In order to guarantee the stability of the coupled system, the solution of the arterial 1D model is initialised until a periodic solution is reached. This is achieved by performing 20 cycles of the circulatory model alone prior to coupling with the heart model, with a flow or pressure profile provided as an inlet boundary condition. For example, a flow or pressure waveform clinically measured at the aortic root can be used as a boundary condition, if available. Otherwise, analytical functions can be prescribed Boileau et al. (2015). For the current study, we considered two different choices: an invasivelyrecorded pressure measurement at the aorta and an analytical function of flow waveform calibrated to match measured peak flow, opening, maximum and closing pressure at the aorta. The coupled system is solved using a fully converging Newton method with maximal number of steps and an absolute norm error reduction of the residual of . To illustrate the performance of the coupling method we consider a system consisting of a 3D EM model of the LV and for the 1D arterial model we analyse different geometries, either consisting of a single arterial segment or composed by more complex networks including the major vessels. To reproduce a realistic case, the LV computational domain is adapted from a patientspecific 3DwholeheartMRI scan collected and postprocessed in a previous study, see Figure 1.
The biomechanical parameters of the solid and fluid models under baseline conditions have been calibrated to match a patientspecific set of in vivo measurements at specific instants of the cardiac cycle for the same subject, see Marx et al. (2020) for further information. The same pressurevolume data is used as for the initialisation of the 3D EM model. All parameters used for the baseline case are summarised in Table 1.
5 Results
5.1 Stability assessment of the coupling method
To analyse the stability of the proposed approach, we compared the results of simulations with a varying time resolution for the 1D blood flow scheme or the 3D LV model, respectively. Taking into account the intrinsic timestep limitation necessary to ensure stability of the 1D blood flow scheme, associated with a CFL condition Vos et al. (2010), we did not consider a time step higher than . In more detail, in the first test case we considered , while keeping the time resolution for the 3D LV model constant at . The circulatory system is represented in this case by one vessel segment of and an RCR Windkessel model for terminal boundary conditions. We emphasise that the coupling is always performed at the time resolution of the 3D EM model. Figure 2 demonstrates that solutions remain unchanged when considering a different time discretisation of the circulatory model. Therefore, our results indicate that we consider the lowest temporal resolution ensuring stability of the 1D numerical scheme.
As a second example, we performed two simulations considering a varying time resolution for the 3D heart model, , keeping the same temporal resolution for the 1D circulatory model. Figure 3 shows that, also in this case, solutions remain unchanged when considering a different time discretisation of the 3D heart model. Thus, we will consider in what follows. We refer to Augustin et al. (2016) for a more detailed convergence testing of the cardiac model alone.
5.2 Prediction of pulse wave propagation effects on cardiac function
To test the ability of the coupled model to predict the effect of pulse wave propagation on cardiac electromechanics we consider three test cases, ranging from an idealised aortic segment with a constant radius to a pathological condition called aortic coarctation and a network composed of 116 arteries.
5.2.1 Idealised aortic segment with constant radius
The first example illustrates the effect of vessel stiffening on haemodynamics and LV function. In particular, we considered a 3D LV solid model coupled to a vessel segment with constant radius and a length of with lumped terminal boundary conditions. For the baseline case we set the vessel wall stiffness described by the Young modulus to . For two test cases with increased stiffness we set the Young modulus to and , respectively. The 1D blood flow model was initialised using an inflow boundary condition. For the simulations in this test case we fixed . Simulation results of the coupled model are shown in Figure 4, where the pressure trace in the LV and at the inlet of the aortic segment is compared for the different stiffness values. It can be observed that an increase in produced an increase in peak pressure and a substantial change in the shape of the pressure trace. This is due to the fact that the stiffening of the vessel wall affects pulse wave propagation in the vessel itself. In particular, the increase of aortic stiffness is associated with a premature rise and decay of the pressure wave and the wellknown phenomenon of peak pressure augmentation Nichols et al. (1991). The increased stiffness of the aortic vessel has a significant impact on LV function, see the PV loop in Figure 4 and Figure 5. The stiffening of the arterial vessel is associated with an increase in endsystolic volume (ESV) at unchanged enddiastolic volume (EDV), leading to a reduction of stroke volume (SV)
5.2.2 Stenotic aortic segment
As a second example, we considered a pathological configuration of aortic coarctation, that is a congenital defect consisting of a local narrowing of the aortic lumen (most commonly in the aortic arch). To this end, we considered an arterial segment of endowed with a stenosis model Jin and Alastruey (2021), with a 30% stenosis halfway of its length. Terminal boundary conditions are given by an RCR Windkessel model. In order to explore the effect of stiffening of the vessel in this configuration, we considered a baseline case with and two test cases with equal to and , respectively. For the simulations in this test case we set for stability considerations. Also in this example, an increase in aortic stiffness is reflected in an augmentation of peak pressure and a variation in the pressure profile, see Figure 6. The effect of vessel stiffness on LV function is shown in the PV loop of Figure 6 consists in a reduction in SV caused by an increase of the ESV.
For the sake of completeness we show the variation of pressure and flow signals along the distance of the 1D vessel for the case . In particular, we show the model predictions of the 3D1D coupled model in Figure 8.
5.2.3 116 aortic segments
In order to explore the effect of bifurcations on pulse wave reflections and their impact on LV function the 3D EM model of the LV was coupled to a more complex model of the arterial system. To this end, as a third example we coupled the 3D model to a network consisting of 116 arterial vessels. Moreover, we explored two different physiological configurations for the 1D model, corresponding to a 25 yearold subject and a 65 yearold subject. We refer to Charlton et al. (2019) for further detail on the 1D networks considered. Concerning the 3D heart model, we considered the same parameter values as in the previous test cases, listed in Table 1, with some exceptions. In particular, the heart cycle period was set to in both configurations, and the parameters of the active stress law were adapted according to the new heart cycle period. We set in this test case as well. Also in this example, the increase in aortic stiffness and systemic vascular resistance associated with healthy ageing is reflected in an augmentation of the peak pressure and a variation in the pressure profile, see Figure 9. The effect of ageing on LV function is also striking, as shown in Figures 10 and 9. In particular, we can see a reduction in SV caused by an increase of the ESV. For this test case we can show the effect of pulse pressure amplification, i.e. the amplification in pulse pressure signal towards more distal locations, associated with reflections at bifurcation sites. For the sake of illustration, in Figure 11 we consider the change in pressure trace from the aortic root to the right brachial artery. As expected, pulse pressure amplification was more pronounced in the younger subject than in the 65yo subject.
6 Discussion
In this methodological study we have developed a numerical scheme that allows to couple an advanced 3D EM model of cardiac function to a 1D model of the human circulation. First, we have demonstrated the feasibility of the proposed approach by coupling a patientspecific 3D LV model to a 1D model of the aorta and shown the ability of the model to correctly predict known physiological behaviours associated with arterial wall stiffening both in a healthy case and in the pathological case of aortic coarctation. Furthermore, we have shown an application to the coupling of the 3D LV model with a complex network composed by 116 arteries, in order to study the effect of bifurcations on pulse wave reflections and their impact on cardiovascular function.
6.1 Prediction of pulse wave propagation effects on cardiac function
The results in Section 5 show an augmentation of peak pressure in the LV and pulse pressure in the aorta and a
substantial variation in the shape of the pressure trace as a consequence of
stiffer vessel walls.
Peak pressure in the LV increases with increasing E, since there is a higher afterload.
In the aorta, the pulse pressure increases with increasing Young’s modulus E, since the aortic model becomes less compliant.
Physiologically, the pulse pressure trace is composed of two additive components: first,
a forward propagating wave generated by LV contraction and second,
a backward propagating wave from the periphery.
From this, the behaviour observed in the simulations is to be expected as the
increased aortic stiffness generates a faster pulse wave propagation.
Hence, this causes a premature arrival of the forward propagating wave and, in turn, of the backward propagating wave,
entailing peak pressure augmentation Nichols et al. (1991).
In addition, the stiffening of the arterial vessel is related to
an increase in ESV.
The main reason for this effect is that a lower compliance
of the vessel, i.e., its ability to distend and increase volume with increasing
transmural pressure, is responsible for a reduction of blood ejected by the
heart chamber, thus producing a greater ESV.
Since EDV is influenced by the cardiac passive behaviour and preload rather than
by afterload, this volume does not vary when the aortic vessel properties change.
Hence, the stiffening of the aortic vessel leads to a reduction in SV.
We have shown that the proposed 3D1D model is capable to reproduce the impact
of increased stiffness and pulsatile load in the circulation and their effect
on heart performance.
These results confirm that such models are particularly suited to explore the effects
of changes in pulse wave propagation, which can hardly be captured using
lumped parameter models Alastruey et al. (2012).
Also, effects of these interactions are very complex to study in vivo
due to the need for high fidelity and simultaneous measurements of cardiac function
and circulation at several locations, together with the capability to assess chamber
interactions Mynard and Smolich (2015).
Hence, 3D1D models are of high relevance to investigate
interactions between vascular waves and cardiac chamber function.
Overall, a coupling to 1D models should be preferred over a coupling to
lumped parameter models for clinical applications that are profoundly influenced
by disrupted blood flow propagation, e.g., aortic coarctation pathologies and
pulmonary arterial hypertension.
6.2 Numerical aspects
Computationally, the 3D1D approach is comparable to 3D0D models, as the extra cost of solving the 1D model is essentially negligible compared to the cost of the 3D model. As such, whenever wave transmission effects are to be investigated, the 3D1D model is preferable as it can be used instead of a 3D0D model as in Augustin et al. (2021) without any computational penalties. Further, we emphasise that the implementation of the coupling framework is not constrained to the specific choice of the 1D solver used for the solution of the blood flow equations in the circulatory system. In personalised applications the extra cost of identifying parameters of a distributed 1D system must be factored in though.
6.3 Limitations
An extension of the study is the use of more complex biventricular and
fourchamber EM cardiac models Strocchi et al. (2020a), in order to consider the interactions among the
heart chambers and thus provide more physiological simulations of the cardiac function.
Also, it is well known that cardiac preload is mainly affected by venous return.
As a consequence, in order to study more accurately the feedback of circulation (including venous return) on the cardiac function, the use of a 1D global model of the human circulation is more appropriate.
Thus, another future perspective consists in coupling the 3D EM cardiac model with closedloop global 1D blood flow models, e.g. Müller and Toro (2014), in order to explore the effect of circulation on preload as well.
We emphasise that the coupling strategy holds for a general framework consisting of multiple cardiac cavities and is prone to different discretisations and geometries for the circulatory model, therefore such extensions are feasible.
Third, this study did not specifically focus on the parameterisation of the
coupled model to patientspecific data. However, the personalisation of
cardiovascular models is of crucial importance for their predictive role,
since the circulatory system has a complex network structure and shows a
significant interindividual variability.
Owing to their distributed nature, 1D models of blood flow involve a larger
amount of parameters that must be calibrated than lumped parameter models, related to geometric aspects,
elastic properties and peripheral impedances of the
network Reymond et al. (2011); Blanco et al. (2015).
This poses serious limitations on the applicability of the use of detailed 1D
circulatory models in clinical applications.
Nonetheless, there are numerous contributions in the scientific literature to cope
with this problem, in terms of simplification of the geometry and reduction of the
parameter set, together with the application of optimisation methods and robust
inverse problem strategies.
For example, the topological complexity of the arterial tree can be optimised by
effectively reducing the number of arterial segments included in the 1D model,
still preserving the key characteristics of flow and pressure
waveforms Epstein et al. (2015); Fossan et al. (2018).
Common approaches also contemplate the rescaling of distal properties and vessel
geometry based on allometric scales and global descriptors related to the geometry and properties of the
network that are easy to measure, such as pulse wave velocity, body size,
and biological age Quarteroni et al. (2017).
The level of detail of the network representing the circulatory system should be carefully considered taking into account the specific clinical application under study.
In addition, recent parameter estimation methods based on data assimilation techniques, i.e., algorithms combining mathematical models with available measurements to improve the accuracy of model predictions and estimate patientspecific parameters, have provided promising results in cardiovascular applications
Lombardi (2014); Caiazzo et al. (2017); Müller et al. (2018).7 Conclusion
We developed a stable strategy to perform a coupling between a 3D EM model of the heart and a 1D blood flow model of the arterial system, based on the resolution of a saddlepoint problem for the volume and pressure in the cavity. We built a personalised coupled 3D1D model of LV and arterial system and showed robustness and accuracy of our scheme in a numerical benchmark. We demonstrated the ability of the coupled system to efficiently simulate physiological response to alterations in features of the circulatory system influencing pulse wave transmission, including aortic stiffening, aortic stenosis and bifurcations causing wave reflections. The additional computational costs associated with the use of 1D models instead of standard 0D models in the coupled system are negligible. As a consequence, the use of our coupled 3D1D model is beneficial for a broad spectrum of clinical applications where wave transmission effects are under study, such as aortic coarctation and pulmonary arterial hypertension.
Acknowledgements
The study received support from the Austrian Science Fund (FWF) grant I2760B30 and from BioTechMed Graz grant ILEARNHEART. Simulations for this study were performed on the Vienna Scientific Cluster (VSC4), which is maintained by the VSC Research Center in collaboration with the Information Technology Solutions of TU Wien. The authors acknowledge the financial support by the University of Graz. The authors acknowledge Dr. Laura Marx, PhD (Medical University of Graz, Austria) for the technical support with the calibration of the 3D EM cardiac model.
References
 Augustin et al. (2016) C. M. Augustin, A. Neic, M. Liebmann, A. J. Prassl, S. A. Niederer, G. Haase, G. Plank, Anatomically accurate high resolution modeling of human whole heart electromechanics: a strongly scalable algebraic multigrid solver method for nonlinear deformation, Journal of Computational Physics 305 (2016) 622–646.
 Niederer et al. (2019) S. A. Niederer, J. Lumens, N. A. Trayanova, Computational models in cardiology, Nature Reviews Cardiology 16 (2) (2019) 100–111.
 DiazZuccarini et al. (2013) V. DiazZuccarini, M. Viceconti, V. Stroetmann, D. Kalra, Roadmap for the digital patient, The Digital Patient Community and the Discipulus Consortium: EU Project Discipulus .
 Chabiniok et al. (2016) R. Chabiniok, V. Y. Wang, M. Hadjicharalambous, L. Asner, J. Lee, M. Sermesant, E. Kuhl, A. A. Young, P. Moireau, M. P. Nash, et al., Multiphysics and multiscale modelling, data–model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics, Interface focus 6 (2) (2016) 20150083.
 CorralAcero et al. (2020) J. CorralAcero, F. Margara, M. Marciniak, C. Rodero, F. Loncaric, Y. Feng, A. Gilbert, J. F. Fernandes, H. A. Bukhari, A. Wajdan, et al., The Digital Twin to enable the vision of precision cardiology, European heart journal 41 (48) (2020) 4556–4564.
 Westerhof et al. (2009) N. Westerhof, J.W. Lankhaar, B. E. Westerhof, The arterial windkessel, Medical & biological engineering & computing 47 (2) (2009) 131–141.
 Segers et al. (2003) P. Segers, N. Stergiopulos, N. Westerhof, P. Wouters, P. Kolh, P. Verdonck, Systemic and pulmonary hemodynamics assessed with a lumpedparameter heartarterial interaction model, J Eng Math 47 (34) (2003) 185–199.
 Segers et al. (2008) P. Segers, E. Rietzschel, M. De Buyzere, N. Stergiopulos, N. Westerhof, L. Van Bortel, T. Gillebert, P. Verdonck, Threeand fourelement Windkessel models: assessment of their fitting performance in a large cohort of healthy middleaged individuals, Proc Inst Mech Eng H 222 (4) (2008) 417–428.
 Formaggia et al. (2003) L. Formaggia, D. Lamponi, A. Quarteroni, Onedimensional models for blood flow in arteries, Journal of Engineering Mathematics 47 (2003) 251–276.
 Reymond et al. (2011) P. Reymond, Y. Bohraus, F. Perren, F. Lazeyras, N. Stergiopulos, Validation of a patientspecific onedimensional model of the systemic arterial tree, American Journal of Physiology  Heart and Circulatory Physiology 301 (3) (2011) H1173–H1182, ISSN 03636135.
 Manganotti et al. (2021) J. Manganotti, F. Caforio, F. Kimmig, P. Moireau, S. Imperiale, Coupling reducedorder blood flow and cardiac models through energyconsistent strategies: modeling and discretization, Advanced Modeling and Simulation in Engineering Sciences 8 (1) (2021) 21.
 Van de Vosse and Stergiopulos (2011) F. N. Van de Vosse, N. Stergiopulos, Pulse wave propagation in the arterial tree, Annual Review of Fluid Mechanics 43 (2011) 467–499.
 SainteMarie et al. (2006) J. SainteMarie, D. Chapelle, R. Cimrman, M. Sorine, Modeling and estimation of the cardiac electromechanical activity, Computers & structures 84 (28) (2006) 1743–1759.
 Kerckhoffs et al. (2007) R. C. Kerckhoffs, M. L. Neal, Q. Gu, J. B. Bassingthwaighte, J. H. Omens, A. D. McCulloch, Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation, Annals of biomedical engineering 35 (1) (2007) 1–18.
 Nordsletten et al. (2011) D. Nordsletten, S. Niederer, M. Nash, P. Hunter, N. Smith, Coupling multiphysics models to cardiac mechanics, Progress in biophysics and molecular biology 104 (13) (2011) 77–88.
 Hirschvogel et al. (2017) M. Hirschvogel, M. Bassilious, L. Jagschies, S. M. Wildhirt, M. W. Gee, A monolithic 3D0D coupled closedloop model of the heart and the vascular system: Experimentbased parameter estimation for patientspecific cardiac mechanics, International journal for numerical methods in biomedical engineering 33 (8) (2017) e2842.
 Regazzoni et al. (2020) F. Regazzoni, M. Salvador, P. C. Africa, M. Fedele, L. Dede, A. Quarteroni, A cardiac electromechanics model coupled with a lumped parameters model for closedloop blood circulation. Part I: model derivation, arXiv preprint arXiv:2011.15040 .
 Stevens et al. (2012) G. R. Stevens, A. GarciaAlvarez, S. Sahni, M. J. Garcia, V. Fuster, J. Sanz, RV dysfunction in pulmonary hypertension is independently related to pulmonary artery stiffness, JACC: Cardiovascular Imaging 5 (4) (2012) 378–387.
 Gurev et al. (2015) V. Gurev, P. Pathmanathan, J.L. Fattebert, H.F. Wen, J. Magerlein, R. A. Gray, D. F. Richards, J. J. Rice, A highresolution computational model of the deforming human heart, Biomechanics and modeling in mechanobiology 14 (4) (2015) 829–849.
 Alastruey et al. (2012) J. Alastruey, K. H. Parker, S. J. Sherwin, Arterial pulse wave haemodynamics, BHR Group  11th International Conferences on Pressure Surges (2012) 401–442.
 Guccione et al. (1995) J. M. Guccione, K. D. Costa, A. D. McCulloch, Finite element stress analysis of left ventricular mechanics in the beating dog heart, Journal of biomechanics 28 (10) (1995) 1167–1177.
 Genet et al. (2014) M. Genet, L. C. Lee, R. Nguyen, H. Haraldsson, G. AcevedoBolton, Z. Zhang, L. Ge, K. Ordovas, S. Kozerke, J. M. Guccione, Distribution of normal human left ventricular myofiber stress at end diastole and end systole: a target for in silico design of heart failure treatments, Journal of Applied Physiology 117 (2) (2014) 142–152.
 Niederer et al. (2011) S. A. Niederer, G. Plank, P. Chinchapatnam, M. Ginks, P. Lamata, K. S. Rhode, C. A. Rinaldi, R. Razavi, N. P. Smith, Lengthdependent tension in the failing heart and the efficacy of cardiac resynchronization therapy, Cardiovascular Research 89 (2) (2011) 336.
 Neic et al. (2017) A. Neic, F. O. Campos, A. J. Prassl, S. A. Niederer, M. J. Bishop, E. J. Vigmond, G. Plank, Efficient computation of electrograms and ECGs in human whole heart simulations using a reactioneikonal model, Journal of Computational Physics 346 (2017) 191–211.
 Vigmond et al. (2008) E. J. Vigmond, R. Weber dos Santos, A. J. Prassl, M. Deo, G. Plank, Solvers for the cardiac bidomain equations., Prog Biophys Mol Biol 96 (13) (2008) 3–18.
 Karabelas et al. (2018) E. Karabelas, M. A. Gsell, C. M. Augustin, L. Marx, A. Neic, A. J. Prassl, L. Goubergrits, T. Kuehne, G. Plank, Towards a computational framework for modeling the impact of aortic coarctations upon left ventricular load, Frontiers in physiology 9 (2018) 538.
 Fernández et al. (2007) M. A. Fernández, J.F. Gerbeau, C. Grandmont, A projection semiimplicit scheme for the coupling of an elastic structure with an incompressible fluid, International Journal for Numerical Methods in Engineering 69 (4) (2007) 794–821.
 Gerbeau (2012) J. Gerbeau, P. Moireau, N. Xiao, M. Astorino, CA Figueroa, D. Chapelle, CA Taylor, Biomech Model Mechanobiol 11 (2012) 1–18.
 Mynard et al. (2012) J. P. Mynard, M. R. Davidson, D. J. Penny, J. J. Smolich, A simple, versatile valve model for use in lumped parameter and onedimensional cardiovascular models, International Journal for Numerical Methods in Biomedical Engineering 28 (67) (2012) 626–641.
 Augustin et al. (2021) C. M. Augustin, M. A. Gsell, E. Karabelas, E. Willemen, F. W. Prinzen, J. Lumens, E. J. Vigmond, G. Plank, A computationally efficient physiologically comprehensive 3D–0D closedloop model of the heart and circulation, Computer Methods in Applied Mechanics and Engineering 386 (2021) 114092.
 Gurev et al. (2011) V. Gurev, T. Lee, J. Constantino, H. Arevalo, N. A. Trayanova, Models of cardiac electromechanics based on individual hearts imaging data, Biomechanics and modeling in mechanobiology 10 (3) (2011) 295–306.
 Balay et al. (2018) S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, a. Dave A. M, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc Users Manual, Tech. Rep. ANL95/11  Revision 3.10, Argonne National Laboratory, 2018.
 Henson and Yang (2002) V. E. Henson, U. M. Yang, BoomerAMG: A parallel algebraic multigrid solver and preconditioner, in: Applied Numerical Mathematics, 2002.
 Boileau et al. (2015) E. Boileau, P. Nithiarasu, P. J. Blanco, L. O. Müller, F. E. Fossan, L. R. Hellevik, W. P. Donders, W. Huberts, M. Willemet, J. Alastruey, A benchmark study of numerical schemes for onedimensional arterial blood flow modelling, International journal for numerical methods in biomedical engineering 31 (10) (2015) e02732.
 Marx et al. (2020) L. Marx, M. A. Gsell, A. Rund, F. Caforio, A. J. Prassl, G. TothGayor, T. Kuehne, C. M. Augustin, G. Plank, Personalization of electromechanical models of the pressureoverloaded left ventricle: fitting of Windkesseltype afterload models, Philosophical Transactions of the Royal Society A 378 (2173) (2020) 20190342.
 Vos et al. (2010) P. E. Vos, S. J. Sherwin, R. M. Kirby, From h to p efficiently: Implementing finite and spectral/hp element methods to achieve optimal performance for lowand highorder discretisations, Journal of Computational Physics 229 (13) (2010) 5161–5181.
 Nichols et al. (1991) W. Nichols, M. O’Rourke, W. L. Kenney, McDonald’s Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles, ed. 3, 1991.
 Jin and Alastruey (2021) W. Jin, J. Alastruey, Arterial pulse wave propagation across stenoses and aneurysms: assessment of onedimensional simulations against threedimensional simulations and in vitro measurements, Journal of the Royal Society Interface 18 (177) (2021) 20200881.
 Charlton et al. (2019) P. H. Charlton, J. Mariscal Harana, S. Vennin, Y. Li, P. Chowienczyk, J. Alastruey, Modeling arterial pulse waves in healthy aging: a database for in silico evaluation of hemodynamics and pulse wave indexes, American Journal of PhysiologyHeart and Circulatory Physiology 317 (5) (2019) H1062–H1085.
 Mynard and Smolich (2015) J. P. Mynard, J. J. Smolich, OneDimensional Haemodynamic Modeling and Wave Dynamics in the Entire Adult Circulation, Annals of Biomedical Engineering (2015) 1–18.
 Strocchi et al. (2020a) M. Strocchi, C. M. Augustin, M. A. Gsell, E. Karabelas, A. Neic, K. Gillette, O. Razeghi, A. J. Prassl, E. J. Vigmond, J. M. Behar, et al., A publicly available virtual cohort of fourchamber heart meshes for cardiac electromechanics simulations, PloS one 15 (6) (2020a) e0235145.
 Müller and Toro (2014) L. O. Müller, E. F. Toro, A global multiscale mathematical model for the human circulation with emphasis on the venous system, International journal for numerical methods in biomedical engineering 30 (7) (2014) 681–725.
 Blanco et al. (2015) P. J. Blanco, S. M. Watanabe, M. R. F. Passos, P. A. Lemos, R. A. Feijóo, An anatomically detailed arterial network model for onedimensional computational hemodynamics, IEEE Transactions on Biomedical Engineering 62 (2) (2015) 736–753.
 Epstein et al. (2015) S. Epstein, M. Willemet, P. J. Chowienczyk, J. Alastruey, Reducing the number of parameters in 1D arterial blood flow modeling: less is more for patientspecific simulations, American Journal of PhysiologyHeart and Circulatory Physiology 309 (1) (2015) H222–H234.
 Fossan et al. (2018) F. E. Fossan, J. MariscalHarana, J. Alastruey, L. R. Hellevik, Optimization of topological complexity for onedimensional arterial blood flow models, Journal of the Royal Society Interface 15 (149) (2018) 20180546.
 Quarteroni et al. (2017) A. Quarteroni, T. Lassila, S. Rossi, R. RuizBaier, Integrated Heart–Coupling multiscale and multiphysics models for the simulation of the cardiac function, Computer Methods in Applied Mechanics and Engineering 314 (2017) 345–407.
 Lombardi (2014) D. Lombardi, Inverse problems in 1D hemodynamics on systemic networks: a sequential approach, Int. J. Numer. Method Biomed. Eng. 30 (2) (2014) 160–179.

Caiazzo et al. (2017)
A. Caiazzo, F. Caforio, G. Montecinos, L. O. Muller, P. J. Blanco, E. F. Toro, Assessment of reducedorder unscented Kalman filter for parameter identification in 1dimensional blood flow models using experimental data, International journal for numerical methods in biomedical engineering 33 (8) (2017) e2843.

Müller et al. (2018)
L. O. Müller, A. Caiazzo, P. J. Blanco, Reducedorder unscented Kalman filter with observations in the frequency domain: application to computational hemodynamics, IEEE Transactions on Biomedical Engineering 66 (5) (2018) 1269–1276.
 ten Tusscher et al. (2004) K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, A. V. Panfilov, A model for human ventricular tissue., American Journal of Physiology. Heart and Circulatory Physiology 286 (2004) H1573–H1589.
 Land and Niederer (2018) S. Land, S. A. Niederer, Influence of atrial contraction dynamics on cardiac function, International journal for numerical methods in biomedical engineering 34 (3) (2018) e2931.
 Strocchi et al. (2020b) M. Strocchi, M. A. Gsell, C. M. Augustin, O. Razeghi, C. H. Roney, A. J. Prassl, E. J. Vigmond, J. M. Behar, J. S. Gould, C. A. Rinaldi, M. J. Bishop, G. Plank, S. A. Niederer, Simulating ventricular systolic motion in a fourchamber heart model with spatially varying robin boundary conditions to model the effect of the pericardium, Journal of Biomechanics 101 (2020b) 109645.
 Klotz et al. (2007) S. Klotz, M. L. Dickstein, D. Burkhoff, A computational method of prediction of the enddiastolic pressurevolume relationship by single beat, Nat Protoc 2 (9) (2007) 2152–8.
 Walker et al. (2005) J. C. Walker, M. B. Ratcliffe, P. Zhang, A. W. Wallace, B. Fata, E. W. Hsu, D. Saloner, J. M. Guccione, MRIbased finiteelement analysis of left ventricular aneurysm, American Journal of PhysiologyHeart and Circulatory Physiology 289 (2) (2005) H692–H700.
 Formaggia et al. (2001) L. Formaggia, J.F. Gerbeau, F. Nobile, A. Quarteroni, On the Coupling of 3D and 1D NavierStokes equations for Flow Problems in Compliant Vessels, Comp. Meth. Appl. Mech. Engrg. 191 (67) (2001) 561–582.
 Westerhof et al. (1971) N. Westerhof, G. Elzinga, P. Sipkema, An artificial arterial system for pumping hearts., Journal of applied physiology 31 (5) (1971) 776–781.
 Alastruey et al. (2008) J. Alastruey, K. H. Parker, J. Peiró, S. J. Sherwin, Lumped Parameter Outflow Models for 1D Blood Flow Simulations: Effect on Pulse Waves and Parameter Estimation, Communications in Computational Physics 4 (2008) 317–336.
 Blanco et al. (2014) P. J. Blanco, S. M. Watanabe, E. A. Dari, M. R. F. Passos, R. A. Feijóo, Blood flow distribution in an anatomically detailed arterial network model: criteria and algorithms, Biomech. Model Mechanobiol. 13 (6) (2014) 1303–1330.
 Votta et al. (2013) E. Votta, T. B. Le, M. Stevanella, L. Fusini, E. G. Caiani, A. Redaelli, F. Sotiropoulos, Toward patientspecific simulations of cardiac valves: stateoftheart and future directions, Journal of biomechanics 46 (2) (2013) 217–228.
 Holzapfel (2000) G. A. Holzapfel, Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons Ltd, Chichester, 2000.
 Deuflhard (2011) P. Deuflhard, Newton methods for nonlinear problems: affine invariance and adaptive algorithms, vol. 35, Springer Science & Business Media, 2011.
 Press et al. (2007) W. H. Press, S. A. Teukolsky, W. T. Vettering, B. P. Flannery, NUMERICAL RECIPES The Art of Scientific Computing Third Edition, 2007.
 Rumpel and Schweizerhof (2003) T. Rumpel, K. Schweizerhof, Volumedependent pressure loading and its influence on the stability of structures, International Journal for Numerical Methods in Engineering 56 (2) (2003) 211–238.
Appendix A Electromechanical PDE model
This section is devoted to the description of the mathematical models considered to describe the most fundamental aspects of cardiovascular function, comprising electrophysiology, passive and active mechanics and haemodynamics.
a.1 Electrophysiology
A recently developed RE model Neic et al. (2017) is considered for the generation of electrical activation sequences serving as a trigger for active stress generation in cardiac tissue. This hybrid RE model combines a standard RD model based on the monodomain equation with an eikonal model. In more detail, we consider the eikonal equation given as
(6) 
where is the gradient with respect to the enddiastolic reference configuration , is a positive function that describes the wavefront arrival time at location , and denote the initial activations at locations . The symmetric positive definite tensor contains the squared velocities associated with the tissue’s eigenaxes , , and . Then, the arrival time function is used in a modified monodomain RD model given as
(7) 
where an arrival time dependent foot current, , is added, in order to mimic subthreshold electrotonic currents and produce a physiological foot of the action potential. The key advantage of this RE model is that it enables to compute activation sequences at much coarser spatial resolutions without being afflicted by the spatial undersampling artefacts leading to conduction slowing or even numerical conduction blocks that are observed in standard RD models. The tenTusscher–Noble–Noble–Panfilov model of the human ventricular myocyte ten Tusscher et al. (2004) is employed to model Ventricular EP.
Parameter  Value  Unit  Description 
Passive biomechanics  
tissue density  
650  bulk modulus  
0.8  stiffness scaling  
5.0  []  fibre strain scaling  
6.0  []  crossfibre inplain strain scaling  
3.0  []  radial strain scaling  
10.0  []  shear strain in fibresheet plane scaling  
2.0  []  shear strain in fibreradial plane scaling  
2.0  []  shear strain in transverse plane scaling  
Active biomechanics  
0.7  onset of contraction  
60.0  membrane potential threshold  
15.0  EM delay  
60  peak isometric tension  
575.0  duration of active contraction  
105.0  baseline time constant of contraction  
90.0  time constant of relaxation  
35.0  []  degree of lengthdependence  
100.0  lengthdependence of upstroke time  
Electrophysiology  
cycle time ()  
AA delay  20.0  interatrial conduction delay  
AV delay  100.0  atrioventricular conduction delay  
VV delay  0.0  interventricular conduction delay  
(0.6, 0.4, 0.2)  conduction velocities  
(0.44, 0.54, 0.54)  conductivities in LV  
1/1400  membrane surfacetovolume ratio  
1  membrane capacitance 
a.2 Passive mechanics of the heart
The deformation of the heart is governed by imposed external loads such as pressure in the cavities or from surrounding tissue and active stresses intrinsically generated during contraction. The cardiac tissue is characterised as a hyperelastic, nearly incompressible, anisotropic material with a nonlinear stressstrain relationship. Mechanical deformation is described by Cauchy’s equation of motion leading to a boundary value problem
(8) 
for , where is the density in the Lagrange reference configuration, is the unknown nodal displacement, is the nodal velocity, is the nodal acceleration, is the deformation gradient, is the second Piola–Kirchhoff stress tensor, represents the body forces and denotes the divergence operator in the reference configuration. We consider as initial conditions
and we set for the sake of simplicity. The boundary of the left ventricular model is decomposed in three components, i.e., , where denotes the endocardium, the epicardium, and the base of the ventricle. At the endocardium normal stress boundary conditions are imposed:
(9) 
with the outer normal vector.
Omnidirectional spring type boundary conditions
constrain the ventricle at the basal cut plane Land and Niederer (2018),
and spatially varying normal Robin boundary conditions
are applied at the epicardium Strocchi et al. (2020b) to simulate the pericardial mechanical constrains.
We consider the following decomposition of the total stress :
(10) 
where and denote the passive and active stresses, respectively. Passive stresses are modelled according to the constitutive law
(11) 
with an invariantbased strainenergy function used to model the anisotropic behaviour of cardiac tissue. Numerous constitutive models have been developed in the literature, ranging from simpler transverselyisotropic models to more complex orthotropic laws. In this article, we consider a hyperelastic and transversely isotropic strainenergy function
(12) 
proposed by Guccione et al. Guccione et al. (1995). Here, the term in the exponent is
(13) 
and is the modified isochoric Green–Lagrange strain tensor, where with . Default values of , , and are used. The parameter is used to fit the LV model to an empirical Klotz relation Klotz et al. (2007) by a combined unloading and reinflation procedure. The bulk modulus , which serves as a penalty parameter to enforce nearly incompressible material behaviour, is set as .
Mechanics boundary conditions at the LV
The different BCs applied to the LV models are summarised in Fig. 12. The springs attached to the aortic rim and at the pericardium are shown in Fig. 12A, as well as the pressure BC in the cavity at the endocardial surface. The pericardial springs penalise displacement only in normal direction and are gradually scaled from the apex to the base. Therefore, the distance in apicobasal direction is used to create a penalty map, see Fig. 12A. To avoid nonphysiological rotation, further springs are attached to the septum, see Fig. 12B. The location of the septal springs is selected automatically by constructing a local coordinate system spanned by the centres of the apical region, the MV and the AV.
a.3 Active mechanics
Following the approach proposed by Genet et al. (2014); Walker et al. (2005), we assume that stresses due to active contraction are orthotropic, with full contractile force along the myocyte fibre orientation and contractile force along the sheet orientation . In more detail, we define the active stress tensor as
(14) 
with the scalar active stress that describes the contractile force. Active stress generation is modelled using a simplified phenomenological contractile model Niederer et al. (2011). Owing to its small number of parameters and its direct relation to clinically measurable quantities such as peak pressure and the maximum rate of rise of pressure, this model is easy to fit and thus suitable for being used in clinical EM modelling studies. In particular, the active stress transient is defined as
(15) 
where
(16) 
and is the onset of contraction, is a nonlinear lengthdependent function in which is the fibre stretch and is the lower limit of fibre stretch below which no further active tension is generated, is the local activation time from Eq. (6), is the EM delay between the onsets of electrical depolarisation and active stress generation, is the peak isometric tension, is the duration of active stress transient, is time constant of contraction, is the baseline time constant of contraction, is the lengthdependence of , is the time constant of relaxation, and is the degree of length dependence. Thus, active stresses are only lengthdependent in this simplified model, but dependence on fibre velocity is neglected.
Appendix B Circulatory system: 1D PDE model
In what follows, the main mathematical aspects of 1D blood flow modelling in the human circulation are detailed.
b.1 Some features of 1D blood flow modelling
The 1D equations of blood flow can be derived by integrating the NavierStokes equations of a generic section of the tube (after assuming axisymmetric flow in a cylindrical tube of radius ), obtaining
(17) 
where is the crosssectional area of the lumen, is the average flux and is the average internal pressure over the crosssection. Note that in the momentum balance a momentum flux correction factor (also known as Coriolis coefficient) is introduced, defined as
An energy inequality that bounds a measure of the energy of the hyperbolic system was derived in Formaggia et al. (2001). Alternatively, it is possible to rewrite the system in terms of the variables . In more detail, if we assume a flat velocity profile, i.e. , we obtain the system
(18) 
where is the crosssectional area of the lumen, is the average axial velocity and is the average internal pressure over the crosssection.
For a constant velocity profile satisfying the noslip condition, the friction force per unit length is , where is a coefficient depending on the blood viscosity and the average axial velocity. As System (18) is composed by two equations and three unknowns (), a closure condition is needed, namely the socalled tube law:
(19) 
where and are parameters depending on the wall stiffness and the wall viscosity, respectively.
A simple, purely elastic tube law consists in defining
with the Young modulus and the wall thickness of the vessel. A more complex tube law is obtained when viscoelastic effects are considered, e.g., by setting
with the wall viscosity. In what follows we consider a viscoelastic tube law, see table 2 for more detail.
In the first test cases proposed in this work we consider an arterial network consisting of a segment of the human upper thoracic aorta, equipped with a 3element Windkessel model as a terminal boundary condition. A more complex network of vessels is also considered in Section 5.2.3, and it is taken from Charlton et al. (2019).
Parameter  Value  Unit  Description 

fluid density  
4e3  blood flow viscosity  
1.1  []  Coriolis coefficient (momentum equation)  
1  vessel wall viscosity 
b.1.1 Boundary conditions
Standard inlet boundary conditions for 1D blood flow models consist of lumped 0D models of the heart or, alternatively, pressure or flow profiles Formaggia et al. (2003); Alastruey et al. (2012); Müller and Toro (2014); Blanco et al. (2015). Every arterial 1D model must be truncated after a certain number of generations of bifurcations, since the relative size of red blood cells to vessel diameter increases, therefore the assumption made in 1D modelling that blood is a continuum and a Newtonian fluid fails. Moreover, contrary to large arteries, fluid resistance dominates over wall compliance and fluid inertia in smaller vessels. For these reasons, linear lumped parameter models are commonly employed to simulate the effect of peripheral vessels on pulse wave propagation in larger arteries Alastruey et al. (2012); Westerhof et al. (1971). Such models are obtained by averaging Eq. (18) over the length of a vessel and considering some simplifying hypotheses, such as neglecting the convective term in the momentum equation. A suitable lumped parameter model for our purposes is composed of one capacitor and two resistors Alastruey et al. (2008). In more detail, the first resistance corresponds to the characteristic aortic impedance. This is connected in series with a parallel combination of the peripheral arterial resistance and compliance . We define as the pressure at which flow in the microcirculation is equal to zero. The resulting model is governed by
(20) 
b.1.2 Numerical implementation
Numerous models have been proposed in recent years for the numerical discretisation of the resulting system of hyperbolic equations, relying on different discretisation approaches mainly based on the finite element method (FEM) Formaggia et al. (2003); Blanco et al. (2014) or finite volume method Müller and Toro (2014). A recent benchmark study on six commonly used numerical schemes for 1D blood flow modelling Boileau et al. (2015) showed a good agreement among these schemes and their capability to capture the main characteristics of flow and pressure waveforms in the large arteries. In this study we used the solver Nektar1D Alastruey et al. (2012), which is based on a FEM discretisation (discontinuous Galerkin) in space and finite difference (explicit secondorder Adams–Bashforth) in time of Equations (18) and (19).
Appendix C Valve dynamics models
The EM model of the heart is coupled to the arterial network via the cardiac valves. Flow through these valves is basically triggered by a positive difference pressure between the two compartments. However, several models can be considered for modelling the dynamics of these valves, ranging from simple diodelike models to more sophisticated approaches, that are able to contemplate pathological conditions such as a stenosis or regurgitation. In the literature there are three basic approaches, ranging from simple diodelike models to more complex valve models Quarteroni et al. (2017), simulating a simplified valve dynamics without explicitly modelling the valve leaflets. However, when it is necessary to capture more complex flow dynamics inside the ventricles that are induced by valve pathologies, it may be required to explicitly model the shape and dynamics of the leaflets. Several approaches have been proposed Votta et al. (2013), but this is still a very open and demanding field of research. The simplest model to simulate the action of a valve corresponds to a diode associated with a resistance (characteristic of the valve), with instant opening and closing depending on the pressure gradient upstream and downstream the valve. This model has a physiological basis, since cardiac valves open and close very rapidly (in a few milliseconds). However, it does not take into account any dynamics on , and it is not possible to model pathological conditions. As a consequence, for this work we considered a more accurate model, namely a lumped parameter model of valve dynamics proposed in Mynard et al. (2012). Its peculiarity consists in a smooth opening and closing dynamics of the valves using a simple approach. In particular, this model is derived from the Bernoulli equation, after neglecting Poiseuilletype viscous losses (proportional to ), since they are small:
(21) 
The Bernoulli resistance is associated with the pressure difference related to convective acceleration and dynamic pressure losses caused by the diverging flow field downstream to the vena contracta, and it reads
(22) 
where is the blood density (usually set equal to ) and is an effective crosssectional area. The blood inertance is related to the pressure difference associated with blood acceleration and it is governed by the following equation:
(23) 
with an effective length, and it is here taken equal to the diameter of . In order to consider valve dynamics, is a state variable and depends on an index of valve state (: closed valve, : open valve) such that
(24) 
Note that, in order to avoid division by zero in Equations (22) and (23), the valve state must be always strictly positive in practice. In order to consider pathological conditions like valve regurgitation, stenosis and varying annulus area , the maximum and minimum effective areas can be expressed as
(25) 
In this study, is considered constant. The coefficients and range within , with and in case of a healthy valve, whereas indicates a regurgitant valve and stands for a stenosed valve. Finally, corresponds to an absent valve (unobstructed conduit), whereas stands for an atretic valve.
The rate of valve opening and closing only depends on two factors, i.e., the instantaneous pressure difference across the valve and the current state of the valve. In particular, it is assumed that the valve opens when exceeds a threshold , and that the rate of opening is governed by
(26) 
with a rate coefficient for valve opening (in ). On the other hand, it is assumed that the valve closes when is lower than a threshold , and that the rate of closing is governed by
(27) 
where is a rate coefficient for valve closing (in ). Model parameters are depicted in Table 3.
Parameter  Value  Unit  Description 
Aortic Valve  
2  rate coefficient for valve opening  
0  threshold pressure difference for valve opening  
1.2  rate coefficient for valve closing  
0  threshold pressure difference for valve closing  
1  []  valve stenosis coefficient  
0  []  valve regurgitation coefficient  
Mitral Valve  
0.2  rate coefficient for valve opening  
0  threshold pressure difference for valve opening  
0.2  rate coefficient for valve closing  
0  threshold pressure difference for valve closing  
1  []  valve stenosis coefficient  
0  []  valve regurgitation coefficient 
Appendix D Coupling strategy
Various approaches can be envisaged to perform the coupling between a (3D) cardiac model and a circulatory model. In practice, the problem consists in finding the new state of deformation as a function of the pressure in the cavity applied as a Neumann boundary condition at the cavitary surface. This unknown pressure needs to be appropriately determined. Two configurations need to be captured by the model:

When all valves are closed, the cavity is in an isovolumetric state, i.e., the muscle enclosing the cavity may deform, but the volume has to remain constant. Consequently, during an isovolumetric phase the pressure in the cavity needs to adjust to the variation over time of active stresses, in order to maintain the cavitary volume constant.

When one valve is open or regurging, there is a change in cavitary volume. In this configuration the pressure is regulated by the state of the circulatory system or of a connected cavity. Therefore, needs to be estimated by matching mechanical deformation and state of the system. In fact, pressure in the cardiovascular system depends on flow, which depends in turn by cardiac deformation. As a consequence, the heart and circulatory models are tightly bidirectionally coupled.
From a mathematical perspective, this coupling can be performed in two ways. The simplest approach is to prescribe explicitly. Then, the coupling is performed during the ejection phase by updating flow and flow rate based on the current prediction on the change in the state of deformation associated with the currently predicted pressure . In particular, the pressure boundary condition in each nonlinear solver step is modified within each iteration . Note that the new prediction is then explicitly prescribed as a Neumann boundary condition. Therefore, this explicit and partitioned approach may introduce some inaccuracies during ejection phases and can lead to instabilities during isovolumetric phases. Such instabilities arise from the difficulty of estimating the change in pressure necessary to keep the volume constant. In fact, a knowledge on cavitary elastance is required to know the relation of the cavity at this given point in time, that is not available. Therefore, iterative estimates are needed to gradually inflate or deflate a cavity to its prescribed volume. However, since the elastance properties of the cavities are highly nonlinear, an overestimation may induce oscillations, whereas an underestimation may induce very slow convergence and prohibitively large numbers of Newton iterations. A more sophisticated approach consists in treating as an unknown. Thus, an additional equation is required to close the system, leading to a saddle point formulation and a block system to be solved (monolithic approach).
Isovolumetric Phases
During the isovolumetric phases we can state that
i.e., the volume of each cavity equals an initial volume and does not change during the isovolumetric phase. Hence, in system (4) the matrix and .
Appendix E Finite Element Formulation
e.1 Variational Formulation
For the sake of simplicity of formulation, we neglect here the acceleration term in Eq. (3) in this section and consider at the stationary counterpart of the boundary value problem (8–9) with (1). Please refer to Augustin et al. (2021) for further detail on the full nonlinear elastodynamics problem. The boundary value problem (8–9) with (1) can be formally expressed in terms of the equations
(28)  
(29) 
which is valid for all vector fields smooth enough and vanishing on the Dirichlet boundary , test functions that are for the cavity and otherwise, with the duality pairing . The first term on the left hand side of the variational Eq. (28) represents the rate of internal mechanical work and is defined as
(30) 
where denotes the second Piola–Kirchhoff stress tensor, see (10), and is the directional derivative of the Green–Lagrange strain tensor, see Augustin et al. (2016); Holzapfel (2000). The second term on the left hand side denotes the weak form of the contribution of pressure loads (28) and is computed using (9)
(31) 
The first term of the coupling equation (29) is obtained from (2) using Nanson’s formula and by
(32) 
The second term of (29) is computed using the 1D blood flow model, see B, for .
e.2 Consistent Linearisation
A Newton–Raphson scheme is employed to solve the nonlinear variational equations (28)–(29), with a finite element approach see Deuflhard (2011). Let us consider a nonlinear and continuously differentiable operator . Then, a solution to can be approximated by
which is evaluated until convergence. In our approach, we assume , , , , and . Thus, we obtain the following linearised saddlepoint problem for each , find such that
(33)  
(34) 
with the updates
(35)  
(36) 
The Gâteaux derivative of (28) with respect to the displacement change update yields the first term in (E.2)
(37) 
and the second term in (E.2), that is given by
(38) 
where we have defined, for the sake of simplicity:
The third term in (E.2) is defined using the Gâteaux derivative of (28) with respect to the pressure change update :
Comments
There are no comments yet.