Coupling reduced-order blood flow and cardiac models through energy-consistent strategies: modeling and discretization

In this work we provide a novel energy-consistent formulation for the classical 1D formulation of blood flow in an arterial segment. The resulting reformulation is shown to be suitable for the coupling with a lumped (0D) model of the heart that incorporates a reduced formulation of the actin-myosin interaction. The coupling being consistent with energy balances, we provide a complete heart-circulation model compatible with thermodynamics hence stable numerically and informative physiologically. These latter two properties are verified by numerical experiments.


Introduction
The importance of reduced-order (RO) models in clinical applications has been extensively assessed in the last years. In particular, RO models are nowadays very widespread in the scientific literature [1][2][3][4][5] concerning cardiovascular applications for patient-specific model predictions. Lumped-parameter zero-dimensional (0D) models-typically Windkessel models [6][7][8]-can provide a general view on the global response, e.g. in pressure and flow, of the whole cardiovascular system or a portion of it [4,9]. Hence, RO models can be used as simplified limit conditions for a more detailed system, e.g. a full threedimensional (3D) heart model [10][11][12][13] and in practice they are often used to represent the circulation upstream and/or downstream of the domain of a higher-order model and to define the relationship between pressure and flow at its boundaries [14][15][16].
RO models used in cardiovascular applications are one-dimensional (1D) models of the blood circulation. The 1D formulation accounts for the effects due to pulse wave transmission and thus enables to predict important markers such as pressure wave velocity. These models have been validated against in vitro [17,18] and in vivo [19,20] measurements and have proven to provide useful insights for the understanding of cardiovascular physiology and pathology. Further, 1D models may be preferred over 0D models when local vascular changes or distributed properties (e.g. tapering, branching, stenoses) are under study, and when the impact of physiological and disrupted wave transmission on the circulation (f.e. the origin and clinical relevance of the dicrotic notch, as in [21]) is investigated. Therefore, 1D models are ideal to simulate blood flow in a single arterial segment [14,20,22,23] or a more structured tree of arteries [16][17][18]24]. In addition, 1D models may also be used to take into account other components (e.g. coronary circulation) [19] and be employed to simulate the global human circulation [25][26][27][28] (in combination with lumped-parameter models for the heart dynamics [29], pulmonary circulation and microvascular beds).
In addition, RO models can also be employed to model the heart dynamics. The simplest approach to model the heart dynamics is based on the use of time-varying elastance heart dynamics models [30]. More accurate RO approaches were also developed, see for instance [31,32]. In the 0D model proposed in [32], the geometry of the left ventricle is considered to be a thick sphere and the dependence on space of the heart motion is only related to the radius of the sphere. In [31], the mechanics of ventricular interaction is based on the assumption of a simplified ventricular composite geometry. In more detail, ventricular geometry is approximated by three thick-walled spherical segments encapsulating the LV and RV cavities.
Due to their reduced computational cost, RO models are well suited for real-time monitoring when coupled with data assimilation strategies [33][34][35]. Furthermore, due to the reduced number of parameters used to describe the reduced dynamics (compared to full 3D models for instance), they are more adequate for the stable solution of inverse problems (IP). However, IP strategies and data assimilation strategies may fail if the forward problem lacks of appropriate stability properties. Typically, the forward model must provide a stable solution, especially when it is solved with various sets of parameters and (noisy) feedback terms. In this regard, energy-preserving schemes are ideal to discretize forward problems, since they are stable with respect to a variation of parameters and they ensure a reliable control on the behavior of the solution [36]. Moreover, energy balance and exchanges turn out to be important physiological markers that should be well approximated and be easily accessed, thus motivating once more the use of RO models that preserve the energy balance intrinsic to the considered modeled phenomena.
Several difficulties must be overcome when deriving energy-preserving RO models and their associated discretization. First at the continuous level, RO models do not necessarily come under a form that is obviously compliant with energy principles (appropriate energy balance may be lost during the model reduction). Hence, such RO models must be adapted or transformed. Then, it is also not obvious to construct a numerical scheme that preserves the continuous energy balance, since advanced RO models are often highly non-linear.
In this work, we extend the previous work of [29] by introducing a coupling strategy involving the heart model of [32] and an arterial segment. On the one hand, we provide a novel energy-consistent formulation for the classical 1D formulation of blood flow in an arterial segment. Of note, although the application envisaged blood circulation, this energy-consistent mathematical framework and the formulation proposed are prone to extension to other non-linear hyperbolic 1D problems, like shallow water equations. The resulting reformulation is shown to be suitable for the coupling with the lumped (0D) model of the heart initially proposed in the work of [32]. On the other hand, the novelty concerning the cardiac model is that we are able to prove the stability both in the continuous and discrete reduced order formulations, hence guaranteeing the consistency with the energy relation described by [36] for a 3D heart model. The coupling being consistent with energy balances, we finally provide a complete heart-circulation model compatible with thermodynamics hence stable numerically and informative physiologically.
The paper is structured as follows: first, in "Partial differential equations of the reduced models" section we present the equations that describe, respectively, the 1D blood flow model, the 0D heart model and the valve dynamics, and we report their energy relation for all the compartments. Then, in "An energy-compliant formulation for the blood flow model" section we detail the changes of variables that are performed in the blood flow model to obtain an energy-compliant formulation and we illustrate its non trivial numerical aspects. The key aspects of the discretization are presented in "Discretization" section. Finally, in "Simulations and results" section we show the results of the numerical simulations and we draw the conclusions and perspectives.

Blood propagation in the aorta
Momentum and continuity equations in their one-dimensional formulation are widely used to model the arterial tree, or a portion of it, and to study pressure and blood flow [19]. The standard formulation that describes blood flow propagation in a vessel is derived from the Navier-Stokes equation by an asymptotic analysis procedure [2,37]. For all time t > 0 we look for a blood flow Q(t, s) and a cross-section A(t, s), along the axis s of the vessel, that are solution of where α vp is a coefficient related to the a-priori assumption on the velocity profile in the vessel, P(A) is the pressure, P ext the exterior pressure, ρ the density of blood (considered as a constant value) and K r the friction parameter. In this work we assume that the velocity profile is flat, therefore The first equation of System (1) represents a reduced form of the continuity equation, whereas the second one corresponds to the momentum conservation and the last one is a relation that accounts for the vessel wall displacement. In particular, it links the change in pressure to the wall deformation and deformation rate, hence the change in the cross-section A. It reads: where β = (4 √ πE h 0 )/3, with E the Young modulus of the vessel and h 0 its thickness, = (2 √ π h 0 ν)/3, with ν viscosity of the wall and A 0 the reference area of the crosssection of the vessel. The analysis of ψ v will be addressed in "Viscosity of the wall" section. Therefore, if not specified, in the following sections ψ corresponds to ψ e . System (1) should be completed with the initial condition A(0, s) = A 0 and Q(0, s) = 0 as well as boundary conditions that are the subject of the forthcoming sections. When and K r are considered equal to zero, System (1) is composed by non-linear hyperbolic equations and discontinuities may appear in time, e.g. shocks, even when smooth data are considered. However, the presence of viscosity smoothens the solution and it is reasonable to assume that its derivative with respect to s exists.
It is possible to show that, for smooth solutions, an energy relation holds. Surprisingly, such relation has not yet been used in the literature to deduce numerical models for the propagation of blood flow, although it is well known that energy-consistent methods provide great benefits in terms of stability, which is a key aspect in the context of multiphysics couplings.
Considering the continuous problem, if we take into account a vessel of length L, then we can define the energy related to (A, Q)-or to (A, u) where u is the blood velocity since Q = Au-as Moreover, we define the instantaneous loss term as Then, a straightforward extension of [37, Lemma 2.2] can be deduced.

Lemma 1 Any smooth solution of System (1) satisfies the conservation property
with P tot defined as

Outflow and inflow conditions
The energy balance (5) shows that the energy defined in (3) is a decreasing function of time if there is no blood flow imposed at the inlet and the outlet. But of course, System (1) should be completed with more realistic boundary conditions, typically relating the input flux (output flux, respectively) or the pressure described in (6) at the inlet (outlet, respectively) with the flux or pressure in other systems. In our case, a simple three-element Windkessel [38] is employed at the outlet, hence we introduce a new unknown, a pressure P c , that satisfies where C c , R per and R c are positive parameters that correspond to a conductance and resistances in the Windkessel terminology. At the inlet, we have where P ar and Q ar stand for the arterial pressure and the arterial flux at the inlet of the aorta, respectively. One can then easily deduce from Lemma 1 an energy balance for the coupled Eqs. (1) and (7). (1,7) satisfies the conservation property

Theorem 2 Any smooth solution of Systems
where In the next two subsections a non-linear reduced-model of the heart and cardiac valves is derived as an efficient parametrized generator of an inlet flow. This reduced model will then be adequately coupled to the one-dimensional model proposed in this work in order to preserve the energy balance.

Heart model
In this section, we describe the chosen reduced heart model proposed in [32]. This model benefits from an appealing mathematical structure while producing accurate pressurevolume loops. It integrates in a system of ODEs a microscopical Huxley-like model of actinmyosin binding with a macroscopical cavity deformation formulation. The dimension reduction relies on a spherical hypothesis and a shell asymptotic derivation. Here, we will prove that an energy balance exists for the reduced model as it is the case for the threedimensional formulation that is derived in [36]. This will allow us to derive a complete energy balance property when the 0D heart model is coupled to the 1D blood flow.

Cardiac mechanics
Following [32], we can derive by an asymptotic procedure a system of ODEs describing a lumped cardiac mechanical model that is geometrically represented by a thick sphere of radius R and thickness d-see Fig. 2. The unknown displacement field is reduced to a radial lumped quantity y, that is a unique scalar variable such that the deformed radius of the sphere R, the thickness in the deformed configuration d and the volume of the deformed cavity V [39], shown in Fig. 2, are given by where R 0 and d 0 are the radius and the thickness of the sphere in the reference configuration, respectively. The system of ODEs also involves variables accounting for the modeling of heart contraction through the active deformation e c and the active stress τ c , which are linked to the global deformation with the rheology pictured in Fig. 2 following the recent formulation proposed in Kimmig et al. [40]. The system is loaded with the ventricular pressure P v . The dynamics reads where ρ 0 is the density in the reference configuration, μ a viscosity parameter, k s a stiffness parameter accounting for passive components of the myosin filament-typically the passive stiffness of the filaments themselves plus the Z-disks and | 0 | is the volume of the myocardium in the reference configuration, which is given by Moreover W p (y) and W v (y,ẏ) are directly inferred from the passive potential and the viscous pseudo-potential of the connective tissue matrix [32]. They are smooth functions that satisfy in particular W p (y) ≥ 0 and W v (y,ẏ)ẏ ≥ 0.
Choosing an isotropic transverse exponential law [32,41] leads to where C 0 and C 2 are some parameters that describe the stiffness of the material and C 1 and C 3 are non-dimensional parameters, while η is the viscosity.
Assuming for now that the active stress τ c is imposed, the dynamics described in System (10) gives a reduced heart dynamical systems of 2-state variables (y, e c ) for which we can derive the heart energy balance. We first introduce the energy that is a combination of the kinetic energy ρ 0 | 0 |ẏ 2 2 , the hyperelastic energy W p and the elastic energy stored in the series element. Then, we introduce the dissipation term (12) and state the following energy relation result.

Theorem 3 Any smooth solution of System (10) satisfies the energy balance
where, at the right-hand side, we have the coupling term with the circulation and the microscopic active stress input.

Microscopic actin-myosin binding model
The chosen model of active contraction balances our need of reasonable complexity with physiological characteristics regarding pressure-volume loops [36,40,42]. In more detail, the microscopic active stress τ c is computed from the first two moments of a Huxley-like formulation of the actin-myosin binding phenomenon. Introducing the active stiffness variable k c , we have the following system ⎧ ⎨ ⎩τ c = −(|ν| + α|ė c |)τ c + n 0 (e c )σ 0 |ν| + + k cėc , where we use the symbol | · | + to denote the positive part. The parameter k 0 denotes the maximum active stiffness parameter and σ 0 the corresponding maximum active stress, α is a time constant, n 0 (e c ) is a function with values in [0, 1] accounting for the Frank-Starling mechanism, and ν(t) = ν([Ca 2+ ](t)) is a function triggering the contraction, typically when [Ca 2+ ] > c th , with c th a given threshold, see Fig. 1.
The state variables (τ c , k c ) of System (14) model the active contraction triggered by the input signal ν. For the energy balance, we introduce the energy stored in our homogenized model of actin-myosin bridges and a dissipative term [36] Combining (13) and System (15), we finally obtain the complete energy balance for the heart model where, from the point of view of the heart model, ν is an input signal and P v an external loading. (15) is deduced by introducing the variable λ c = τ c / √ k c [36]. Then, as an intermediate step it can be shown that λ c satisfies the following ODE:

Remark 5 The energy relation in System
The energy relation is finally obtained by multiplying (16) by | 0 | λ c .

Cardiac valve models and energy relation for the complete system
The inlet and outlet of the ventricular model are represented, respectively, by the atrioventricular and the aortic valve, as shown in Fig. 2. To represent the fact that they may be open or closed, valves are modeled as diodes. Therefore, using Kirchhoff's circuit laws, we get where the first equation is associated with the atrioventricular valve, whereas the second one relates to the aortic valve. In more detail, P v , P at and P ar represent the pressure in the ventricle, in the atria and at the inlet of the aorta, respectively, whereas Q v and Q ar correspond to the ejected blood flow throughout the ventricle and at the inlet of the aorta, respectively. Finally, K ar , K iso and K at represent the resistances of the valves [43], as depicted in Fig. 2 Now, if we multiply the first equation of System (17) by P v , we multiply the second equation by P ar , we sum them and we define we can obtain an energy relation for the cardiac valve formulation.

Theorem 6 Any smooth solution of System (17) satisfies the conservation property
where we have, at the right-hand side, the input pressure term and the coupling term that take into account both the cardiac and the arterial contribution.
Finally, using Theorems 2 to 6, we are able to retrieve the global energy relation that takes into account the contribution of the heart (including the microscopic modeling of the actin-myosin binding), the valves and the arterial segment.

An energy-compliant formulation for the blood flow model
In order to obtain an energy-preserving scheme for the blood flow model, which satisfies a discrete counterpart of Theorem 1, we will introduce three variational formulations. The first one corresponds to the standard formulation that one obtains directly from System (1). This formulation has A and u as principal unknowns, where A is the crosssection of the aorta defined as A = Q/u and u is the blood velocity. Then, we introduce a second formulation that uses as a primary unknown the radius of the aorta R = √ A/ √ π and u. This formulation is straightforwardly deduced from the first one and is a convenient intermediate step, since it introduces several simplifications. From these intermediate changes we deduce the last formulation that is written in the unknowns where ϕ(R) is a smooth bijective function from R + to I ⊂ R that we define later. This change of variables has the main advantage to provide an "energy-compliant" discretization, as it will be shown in the next chapter. As the reader will see, the energy is a quadratic functional of the new variables (Φ, v) .

Variational formulation in (A, u)
As a first step, we substitute Q with Au in (1) and assume that = 0, hence ψ v = 0. After some algebraic manipulations we obtain that (1) is equivalent to the system Note that we have rewritten System (1) in a specific form adapted to the derivation of the energy balance (in fact System (21) is obtained following the proof of Lemma 1 provided in [37]). Indeed, multiplying the second equation of System (21) by u, one can see that These two equalities are, in fact, essential to prove the energy relation of Lemma 1. The objective is now to derive a weak formulation of System (21). Concerning the first equation, we multiply it by a space-dependent test functionΦ and we integrate in space, obtaining where (·, ·) is the L 2 -scalar product in (0, L). We now focus on the second equation of System (21) and we repeat the procedure performed above, multiplying each term by a space-dependent test functionṽ. After some manipulations, we obtain where a(u; ·, ·) is bilinear in its two last arguments but non-linear in u and is given by and the non-linear functional g is defined as Observe that, by construction, g is linear inṽ and includes only boundary terms. Up to this point, the weak formulation of the problem described in System (21) is Finally, observe that if we substitutẽ in System (27), we can easily retrieve the energy relation presented in Lemma 1. Indeed, thanks to (22), one can see that Moreover, we recover the energy relation Here E ar (t) is the total energy of the 1D model and D ar (t) represents the dissipative term with K r ≥ 0.

An intermediate formulation in (R, u)
In order to obtain a formulation that leads to the achievement of the energy preservation at a discrete level, we construct an intermediate form of System (21). This formulation is obtained by replacing A with πR 2 , where R represents the radius of the lumen. The unknowns become u and R. The first equation of System (21) is now described as Then, the first term in (24) can be rewritten, substitutingṽ ←ṽ/R, as Moreover, one can see that Note that the substitutionṽ ←ṽ/R does not lead to any issue, since we consider solutions with R > 0 at any time and position. Finally, collecting the four expressions above, one can obtain a formulation with R and u as primary unknowns. It reads It is worth noticing that the product Ru appears "almost" naturally and it is therefore tempting to define v := Ru as a new variable. It becomes even more obvious that this choice is suitable by looking at the energy density, defined as This is precisely what motivates the introduction of the next formulation. Moreover, it is worth mentioning that now so ψ(πR 2 ) is linear with respect to the unknown R and we will see in the next sections that (πR 2 ) is a third-order polynomial and this will simplify its analysis.

Remark 8
The change of variable A = πR 2 is still meaningful even if the 1D hemodynamic model does not assume a perfect circle for the geometry of the cross-section. What matters here is that the new variable R depends on the square root of A. Obviously, the introduction of the factor π is natural to obtain a physical meaning for this new variable since, in practice, arterial cross-sections are almost circular.

Variational formulation in (Φ, v)
A change of variables has to be made in order to demonstrate that the scheme is energypreserving after time discretization. More precisely, time discretization can easily deal with energies that involve quadratic terms of the unknowns. However, the energy density described in (33) is not a quadratic term of the unknowns (R, u), but we can see that the first contribution is a quadratic term of v := Ru.
Therefore, we propose to use v as a main unknown. A first naive choice is then to set ϕ(R) equal to (πR 2 ), where (πR 2 ) is defined as and set Φ ≡ ϕ(R) as the other main unknown. However, we show in "Variational formulation in (Φ, v)" section that this choice is not convenient, since ϕ(·) would not be a bijective function from R + to R + . Instead, we define Before studying in more detail the impact of the choice described in System (36) (in particular the bijectivity of the function ϕ), we formally give the variational formulation associated with the new couple of unknowns (v, Φ), where Φ := ϕ(R). Assuming for now that ϕ is bijective, we define the reciprocal function r(Φ) := ϕ −1 (Φ). Then, each term of the second equation of System (32) can be modified as follows: i. The term involving the time derivative reads ii. The non-linear transport term reads , v,ṽ (38) whereã is now a trilinear form. Such reformulation will lead to the choice of an adapted space discretization that preserves, for all sufficiently smooth functions v and u, the propertỹ and in particular for v. Finally, the boundary term g is given by where, for simplicity, we assume that P tot is given. Of course, if more general boundary conditions are considered, g must be modified accordingly.
Using all the expressions above, we obtain the following equation (corresponding to the second equation of System (32)) where, for the sake of clarity, we have written r instead of r(Φ). The first term in (42) clearly shows how the introduction of v simplifies the dynamic behavior of the equation and it will help at the discrete level to demonstrate the energy preservation. Now we deal with the first equation of System (32) in which we use as a test functionΦ ← ξ (R)Φ, with We show in "Variational formulation in (Φ, v)" section that this function is smooth and positive. We obtain If we focus on the first term in (44), we can observe that Now observe that, by definition, ϕ(R) = ± (πR 2 ). Thus, the term above can be rewritten as Since by definition we have Φ = ϕ(R) and R = r(Φ), we can write where again we use the convention r ≡ r(Φ). At this point, the formulation reads One can see in System (47) an apparent lack of symmetry. Indeed, one could expect the second term in the first equation to be equal to the third term in the second equation. This is true however, since we have, using (43), This observation is fundamental to obtain an energy estimate. To summarize, we have deduced from the dynamics (32) the following formulation: for all (Φ,ṽ) sufficiently smooth find, for all with the following initial data This is what we call the energy-compliant variational formulation. At the continuous level, the energy is easily obtained by choosingΦ = Φ andṽ = v. This simple choice of test functions to deduce the energy relation at the continuous level will help in achieving the same energy relation property at a discrete level.

Remark 9
The formulation of System (48) can be obtained for other tube laws ψ(A). However, some properties should be satisfied by the function ψ. In particular, ψ must be at least continuous and

Strong formulation
For the sake of completeness, we show the strong formulation of System (48). Choosing a smooth test function with compact support in [0, L], one can show, using integration by parts, that the following partial differential equations hold: Then, choosing a smooth test function in [0, L] vanishing at the boundaries in System (48) and using integration by parts for System (50), one can deduce the following boundary conditions:

Analysis of the function ϕ(R)
In this section we provide further details on the properties of the function ϕ(R). The definition in System (36) is motivated by the expression of = (πR 2 ) that is rewritten below: where-in this section-R 0 = √ A 0 / √ π is the reference radius of the cross-section. The behavior of this function is shown in Fig. 3. It is straightforward to see that this function, as well as its square root, is not bijective. However, using System (36), the function ϕ(R) is then given by In Fig. 3 we can observe the comparison between (πR 2 ) and ϕ(R). For every R ≥ R 0 the two functions coincide, whereas for R < R 0 they are opposite. However, we can also see that ϕ is bijective from R to some interval I satisfying R + ⊂ I ⊂ R. Moreover, it is easy to prove the following Property.

Analysis of the function ξ(R)
We now focus on the property of the function ξ (R) that is defined by ξ (R) = ψ(πR 2 )/ϕ(R).
In particular we want to check whether the function is smooth and bounded. This is not true because ϕ(R) vanishes and, as one can see in Fig. 3 and in (52), this happens at R = R 0 where-in this section-R 0 = √ A 0 / √ π is the reference radius of the cross-section. Using Eqs. (34) and (52) one can compute . We see in Fig. 3 that ξ (R) has no singularity, it is smooth, strictly positive and monotonically decaying. This result is summarized in the following Property.

Extension of the model in a non-physiological range
There is an equivalence between System (27) and System (48). More precisely, we can state the following Theorem.
Although the bound defined in (53) is expected physiologically, after space discretization there is no guarantee that such property holds intrinsically at any time and any point. Therefore, we propose to modify System (48) for a non-physiological range, e.g. close to R 0, or equivalently, Φ Φ min . More precisely, r(Φ) is not defined for Φ taking smaller values than Φ min . To circumvent this problem we introduce, for a given > 0-a relaxation parameter-the function r , defined by where (a, b) ∈ R 2 are only defined by the constraint that r ∈ C 1 (R). In more detail, one needs to check that hence one can compute that and then a = e bΦ r(Φ ).
The main advantage of using r instead of r is that r is a bijective function from R to R + \ {0}. Hence, at the discrete level for any value of the unknown Φ we are able to compute a corresponding aortic radius R. In this process we introduce-mathematically speaking-a modeling error with respect to System (27). Nevertheless, we have the following straightforward result.
Note that we can choose small enough so that the range of values A(x, t) ∈ (0, A ) for which the mathematical equivalence with System (27) is not satisfied can be set as desired.
In particular, considering the application to hemodynamics, this interval can be chosen so that a solution of System (27) with values A < A is outside the validity of the the tube law described in (2).

Viscosity of the wall
In "One-dimensional blood flow model" section we introduced the third equation of System (1) that relates the pressure with the strain and strain rate of the wall. In particular, it takes into account the velocity of radial displacements [34] thanks to the term ψ v that was assumed to vanish in "Variational formulation in (A, u)" section in order to derive the energy-compliant variational formulation. In this section we address the treatment of this term, ψ v , through the change of variables introduced in "An energy-compliant formulation for the blood flow model" section. Starting from System (21), we have Since we have already dealt in the previous section with the first term, related to ψ e , we focus now on the last one of the equation above, related to ψ v . Starting from (56) and using the first equation of System (21), we obtain This motivates the introduction of the non-linear form c(·; ·, ·) defined by Taking into account the manipulations performed in "An intermediate formulation in (R, u)" and "Variational formulation in (Φ, v)" sections, one can show that Then, it can be shown that the second equation of System (48) can be replaced by Note that the boundary terms in (59) are indeed taken into account, sinceg-defined in (41)-involves the total pressure that is given by (6) and now reads

Outflow conditions, inflow conditions and energy relation
In order to complete the weak formulation of the problem given in System (48), the outflow and inflow conditions need to be specified. This is done by expanding the termg using the coupling condition described in (7) at the outlet, whereas at the inlet we use where P ar (t) and Q ar (t) are the arterial pressure and the arterial flow, respectively. We obtain the following system of equations: Note that a similar energy identity to the one given in Theorem 2 can be derived for this system, as we state below.

Theorem 14 Any smooth solution of System (60) satisfies the conservation property
where, and, Note that System (14) can be easily used with or without coupling with the reduced heart model. Hence, we consider two cases: • Case 1: Imposed inlet flux; In this case the arterial pressure P ar is considered as a new unknown, namely a Lagrange multiplier for the constraint π R(0) v(0) = Q ar . • Case 2: Coupling with the reduced heart model. System (60) should then be completed with Eqs. (10), (14) and (17), that describe the reduced-order cardiac mechanics, the microscopic actin-myosin binding model and the valve model, respectively. Note that in this model Q ar is an unknown that can be straightforwardly substituted in (17) using the relation Q ar = π R(0) v(0).

Time scheme for the blood flow model
In order to obtain the time discretization of the scheme, we assume a given sequence of time instants {t n } n∈N such that t n+1 > t n and we define the time step as t n := t n+1 − t n . Moreover, we define the half time sequence as t n+ 1 2 := t n + t n /2. The scheme proposed below is formally an implicit second-order time discretization scheme: we consider System (48), we rewrite it at time t n+1/2 and we approximate all the terms using a centered finite difference, i. e. for every n Moreover, as in [36], we introduce some intermediate unknowns, i.e. R n+ 1 2 and ξ n+ 1 2 , as follows : With these considerations, we obtain the following semi-discrete problem are not completely characterized yet. These quantities appear below when considering the time discretization of the valve model.
Once (62) is obtained, we have to check if the semi-discrete scheme preserves the total energy, as given in "Outflow conditions, inflow conditions and energy relation" section for the continuous domain. This is the purpose of the following Section.

Semi-discrete energy relation
Energy preservation can be proven rather simply by following the strategy performed at the continuous level. We substitute the test functions with the proper variables as explained in "Variational formulation in (Φ, v)" section. More precisely, we set v = v n+ 1 2 We observe that where, by definition u 2 = (u, u) is the L 2 (0, L)-norm. Identically, we have v n+1 − v n t n , v n+ 1 The first two equations of System (62) become Now, we can observe that the second term of the continuity equation and the third term of the momentum equation are equal (up to a change of sign). Thus, the latter can be substituted with ( Φ n+1 2 − Φ n 2 )/ t n , which gives The equation above can be further simplified by noticing that thanks to (39). Then, we can define, and, and prove, by using a telescopic sum and by adequately dealing with the boundary terms, a semi-discrete equivalent form of (61), showing that Theorem 14 has a counterpart at the semi-discrete level.

Theorem 15 Any solution of System (62) satisfies the following conservation property for all n ∈ N,
Note that the energy relation obtained in Theorem 15 holds with a time step that may vary between each iteration. This property is fundamental since, in practice, cardiac models often adapt the time step to deal with the abrupt changes of phase due the opening and closure of the aortic valve.

Time scheme for the heart model
The discretization of the heart model was performed following a similar approach to the one described in the previous section. This allows to obtain, as for the arterial model, a discrete equivalent of the conservation property described in Section . The time discretization method used in [36] is adapted to the zero-dimensional formulation of the cavity with the additional introduction of the so-called Hilber-Hughes-Taylor (HHT) Method scheme proposed in [44] for the treatment of kinetics variables and inertial terms. The discrete velocityẏ n and the displacement y n are related by the use of an auxiliary variableÿ n -that stands for an approximation of the acceleration. We proceed discretizing the equations of System (10) at time t n+ 1 2 except for the ventricular, atrial and arterial pressures that are interpolated using a θ -scheme (P n+θ . Deviations from the classical implicit mid-point scheme used in "Time scheme for the blood flow model" section are introduced to generate controlled numerical dissipation terms that damp out purely numerical high frequencies modes that cannot be properly resolved with a finite time step-see also "Spurious high frequencies filtering" section. The HHT kinematics equations are and the dynamics equations then read where we use the adapted [45] energy-preserving non-linear choice ∂W p ∂y (y n+ 1 2 ) i fy n+1 = y n = y n+ 1 2 , and additionally ∂V ∂y (y n+ 1 2 ) i fy n+1 = y n = y n+ 1 2 .
Note that in practice, the expressions for y n+1 = y n are implemented using a series development of W p and V , respectively, to avoid numerical rounding errors. Further, is computed from the already proposed energy-balanced time discretization of System (14) and (16) Since the function n 0 (·) has maximum value 1, it is straightforward to prove [36], from backward Euler time discretization, the a priori bound 0 < k n c ≤ k 0 for all n ≥ 1 if it is satisfied at n = 0. Moreover, the time discretization of the variable λ n is consistent with the dynamics of λ c = τ c / √ k c given in Remark 5. In order to retrieve a discrete equivalent of the discrete energy balance obtained in Section , we multiply the second equation of System (68) by (y n+1 − y n )/ t and the third one by | 0 | (e n+1 c − e n c )/ t n , then we sum, obtaining where we have introduced the discrete energy the discrete dissipation term with E n c := Note that (71), (72) and (76) are the time-discrete equivalent of (11), (12) and (13)

Time scheme for the valves and the energy relation for the complete semi-discrete system
The discretization strategy applied for the valve formulation is the same that we showed for the arterial and ventricular model. In particular, we use an implicit mid-point rule, i.e.
and Q ar ∼ Q System (17) is then discretized as follows: where the discrete energy for the valve model reads the dissipation is and the numerical dissipation reads

The discrete energy relation for the global system
Finally, summing Eqs. [(64), (77) and (79)], we obtain the semi-discrete conservation property for the global system.

Corollary 17 Any solution of System (62), (68),(69) and (78) satisfies the following conservation property for all n ∈ N,
Corollary 17 shows that without a source term, i.e when ν n+1 and P n+ 1 2 at vanish, the energy-that is a norm for the solution-is decaying. This is the expected stability property of a robust time discretization.

Space discretization of the blood flow model
The space discretization is rather simple and does not represent a main issue for the global formulation. However, there are some terms, in the aortic model formulation, that have to be treated carefully when choosing the space discretization method. This is the reason why, in this section, we will only present the space discretization of the arterial model without dealing with how the total discretization of the aorta model couples with the other elements of the model. Indeed, these couplings are straightforward from what has been already explained.
First, we introduce a finite dimensional subspace of H 1 (0, L) of continuous functions that is denoted by V h and we assume an interpolation operator I h : C 0 ([0, L]) → V h as given. For each n ∈ N, we look for the solutions , v by v h andṽ byṽ h . Moreover, we have introduced two other notations: i. We use the notation (·, ·) h to represent an approximation of the scalar product in L 2 (0, L) by quadrature formulae. In particular (·, ·) h is a positive definite bilinear form and is equivalent to the L 2 (0, L)-norm in V h . In what follows, for any u h ∈ V h we define u h 2 h := (u h , u h ) h . ii. The critical point from the energetic point of view in performing the space discretization of our scheme can be found in the definition of the term a h . As stated in "Time scheme for the blood flow model" section for the energy conservation, this term has to satisfy the property defined in (39). One way to guarantee this condition is to consider the exact-i.e. we compute the integrals exactly-trilinear formã. That is why the form a h is defined by With this definition and thanks to (39), it is possible to check that, for all Even though the introduction of the interpolation operator in the definition of u involved in a h given by (82) seems unnecessary, one can observe from the structure of the trilinear formã that its last argument would not be polynomial without interpolation. This implies that the application of a standard quadrature method (e.g. the Gauss integration) would not give an exact integration property, and therefore the property stated by (83) may be lost. Using (83), it is straightforward to prove that the energy relation stated in Theorem 15 can be extended to the fully-discrete case. We obtain the following result.

Theorem 18 Any solution of System (62) satisfies the following conservation property for
and we use P 1 -finite elements. The space V h is spanned by the Lagrange nodal basis functions {w j } J j=0 that satisfy the property w j (s i ) = δ ij . The interpolation operator then reads The scalar product (·, ·) h is defined using the trapezoidal rule, while the termã(ṽ h , v h , u h ) should be computed exactly to ensure energy preservation. To do so, it is sufficient to use the Simpson quadrature method in each element for the underlying integration, sincethanks to the use of the interpolation operator-the integrand is a second-order polynomial in each element.

Simulations and results
In order to show the results of our work, we consider three different model settings: • The uncoupled aortic model is composed by a single 1D straight vessel with a homogeneous circular cross-section representing the upper thoracic aorta and by a Windkessel RCR model that takes into account the impedance and compliance of all the remaining vessels at the periphery. In this test case the inlet blood flow is imposed, as shown in Fig. 4. This configuration allows us to study the behavior of the arterial model alone. • The uncoupled cardiac model, depicted in Fig. 5, is composed by the reduced cardiac model described in the previous sections and by a lumped parameter model as outlet boundary condition that represents the entire circulation. More precisely, we have reduced the 1D model of the aorta into an RC-model, the parameters of the resulting 0D model being the equivalent resistance and compliance of the 1D vessel combined with the RCR boundary model parameters. The equivalent resistance and compliance of the 1D vessel are obtained by (formulas readapted from [46]) where μ b is the blood viscosity. This setting allows to show the outcome of the cardiovascular model when the circulation is represented by a simple Windkessel model. • The fully coupled model studied in the present work, with the reduced-dimensional model of the heart and the aortic model coupled through a transmission condition which includes the valve as represented in Fig. 2.
The simulation results are divided into two main sections. In "Numerical validation" section we present the numerical validation of the proposed formulation, starting from the comparison between the results obtained using the uncoupled aortic model and the results of other numerical schemes. Then, we verify that the energy relation is indeed satisfied at the numerical level and discuss its beneficial effect on the computations. In "Physiological outcomes of the coupling" section we highlight the importance of the coupling involving two case studies: dicrotic notch and physiological ageing. The first case study is considered to demonstrate the need of having a one dimensional model as an outlet boundary for the heart to capture all the important features of pressure waves in the larger elastic arteries (e.g. the aorta), whereas in the second one we show how the simulations of specific conditions of the circulation (e.g. stiffening of the vessels due to  ageing) benefit from the presence of a cardiac model whose behavior depends on the arterial conditions. Both cases reflect the natural interplay between the components of the cardiovascular system [1] .

Benchmark for the uncoupled aortic model
In order to support the validity of the developed method, our scheme is tested in one of the benchmarks presented in [18]. In this work, different numerical schemes are compared for 1D arterial modeling in a set of test cases and they made the results freely available. In the following example we take into account the third benchmark configuration [18], that corresponds to the upper thoracic aorta (setting presented in Fig. 4). For this particular test case the authors considered a single uniform vessel represented by a one-dimensional model, that relies on the classical blood flow equations described in System 1 in the equivalent (A, u) formulation-the state variables are the cross section and the velocityand an RCR model as outlet boundary condition. More precisely, in this benchmark, Eqs.
(1), (7) and (8) are solved with Q ar (t) imposed. The parameters of the model and the inflow boundary condition used in [18] were taken from [23] and are reported, respectively, in Table 1 and Fig. 4. We considered these same parameters and inlet conditions, as shown in Fig. 6, and we compared our results with those of other numerical schemes available in the datasets of [18]. In Fig. 6 it is possible to observe that the pressure curves obtained by the numerical schemes reported are consistent, hence validates our implementation. This shows that our formulation is consistent with well-known discretizations. At the price of a higher computational costdue to the implicit nature-the proposed scheme has, in addition, the advantage to offer a provably stable numerical coupling with a numerical approximation of a non-linear description of the heart.

Spurious high frequencies filtering
In the discretization section, we have introduced two parameters α and θ associated with a controlled artificial viscosity of order δt in order to damp possible undesired oscillations (when α > 0 and θ > 1/2). We emphasize that these oscillations are not due to numerical instability, but to energy exchanges between model compartments introduced by the time discretization procedure, in particular the conservative part involving mid-points. We denote these oscillations spurious high frequencies. To illustrate this effect, we show in Fig. 7 the time evolution over a cardiac cycle of the ventricular pressure P n v , the aortic valve flux K −1 ar |P n v −P n ar | + and the proximal arterial pressure P n ar that are computed with the parameters chosen as in [47]. We compare three different configurations, (α, θ) = (0, 0.5) (no artificial viscosity), (α, θ ) = (0, 0.75) (artificial viscosity only on P v ) and (α, θ ) = (1, 0.75) (artificial viscosity onẏ and P v ). It can be seen that the latter case avoids all kinds of spurious frequencies while preserving, qualitatively, the consistency of the approximation and the cost of the numerical scheme. Of note, the artificial viscosity decreases the formal order of accuracy of the scheme, an acceptable price to pay as our time-step is already rather small to account for the stiff parts of the pressure variation.

Numerical validation of the energy balance behind the coupling strategy
As an illustration of the validity of the energy-preserving numerical methods developed in this work, we present the numerical evaluation of the energy conservation residuals of the three main model elements (66), (77) and (79). For that, we define the residuals as Their time evolution is presented in Fig. 8 over several heart cycles for the simulation with the parameters inspired from [47] and whose results are presented in Figs. 9 and 10. We can notice that the magnitude of the residuals is 12 orders of magnitude lower than that of the individual energy fluxes that compose them. This result thus validates the fact that the energy balances are indeed satisfied at the discrete level.  Having an energy-preserving formulation for the discretization of the model equations ensures that we can make sense of the individual terms of the energy balance. On the one hand, it allows to get a better picture on the solution of the model equations by analyzing the fluxes of energy between the different elements of the model. We present in Fig. 9 some of the energy fluxes that can be computed. We define the following notations In the numerical results, we see how the contraction of the atrium creates an initial flux of energy P in,at towards the left ventricle. Then, the chemical energy flux P ATP provided at the micro level is transformed into the mechanical powerẆ by the left ventricle. It can be noted that the model predicts an energetic yield t+T t |Ẇ(τ )| + dτ / t+T t P ATP (τ ) dτ of 26.2%, which lies in the physiological range: 25 to 35% [48,49]. The average power produced by the ventricle over a heartbeat is 1.02 W, which is in accordance with the values obtained experimentally [47,50]. Moreover, the peak of mechanical power reaches 10.4 W, which is consistent with the evaluation of this quantity from other mechanical models for patients suffering from aortic valve disease [51]. The mechanical power is then transferred with little dissipation through the valves into the aorta that receives the influx of energy P ar Q ar . On the other hand, being able to compute the individual energy fluxes inside the model is of great interest from an application point of view, since some of the terms of the energy balance are indeed used by medical doctors. For instance, anesthetists are interested in the work developed by the left ventricle for the monitoring of the heart [47] during surgical intervention, wheres cardiologists see this same parameter as a potential biomarker for the evaluation of left ventricle dysfunctions [50] and the myocardial efficiency is considered as a relevant indicator to assess the state of patients having aortic or valve pathologies [52,53].

Physiological outcomes of the coupling
In this section, we examine the physiological interest of coupling the arterial and the cardiac model. The three model settings described above are involved to highlight the differences in the results when one of the two components, i.e. the cardiac model or the one-dimensional vessel, is not taken into account-the fully coupled model allowing to study the effect of the mutual interplay between the heart and the arterial network.

Importance of the downstream circulation for the heart: the dicrotic notch
The physiological arterial pressure curve shows two main parts: the systolic and the diastolic phase. During the first phase, the heart ejects the blood into the aortic root and the arterial pressure increases rapidly and reaches a peak, known as the systolic pressure value. Then, the pressure starts falling but it is interrupted by an incisura, known as the dicrotic notch, that happens at the time of the passage to the diastolic phase and causes a second peak. The pressure then continues its downslope to its minimum, the diastolic pressure [54,55]. Any change in the shape of this curve represents a modification in the vessel condition, that is why it is important to be able to properly reproduce these features. We focus here on the reproduction of the dicrotic notch. For that purpose, we perform a simulation using the uncoupled cardiac model and the fully coupled model and we compare the results.
From the results presented in Fig. 10 we can observe that only the pressure curve obtained with the fully coupled model (black curve in the center) shows a dicrotic notch. Moreover, this was obtained with an inflow condition (in black in the left box), coming from the heart, which does not present a backflow. This confirms the latest findings of [21] and suggests that the dicrotic notch is not caused by the presence of the backflow but can be observed when wave propagation phenomena are represented. In addition, although both settings are able to reproduce a physiological shape for the curves, there are more differences between the two results beyond the reproducibility of the dicrotic notch. First, we can observe that the peak of the blood flow is slightly increased when the lumped-parameter model is used to represent the full circulation. Moreover, from the pressure-volume loops (in the right box) we can notice a distinction in shape between the two ejections. This suggests that the systolic blood pressure reaches a lower peak value when the circulation is completely represented by a lumped-parameter model. The last consideration is confirmed even more clearly in the pressure curves, since the lumpedparameter model is able to fit the diastolic decay of pressure but it is not sufficient to properly reproduce the systolic peak [46] and, more generally, the systolic phase. In fact, only when a distributed aortic model is considered it is possible to obtain a more physiological wave showing a noticeable dicrotic notch. These results suggest that when a cardiovascular model is used to analyze phenomena strongly associated with the ventriculo-arterial interaction, lumped-parameter models are not sufficient to represent the downstream circulation and it is necessary to employ a higher-dimensional model (e.g. a 1D-model) that is able to capture the influence of pulse wave transmission within the circulation.

Importance of the heart for the aorta: ageing
During ageing, the vessel walls undergo a degeneration of elastin fibers, a decrease in smooth muscle and an increase in collagen. These changes cause the stiffening of the arteries and in particular of the aorta. As a result, it is observed that the systolic peak pressure increases and the aortic and left ventricular late systolic pressure augment, whereas the aortic blood flow peak and the diastolic pressure decrease. These features are indeed commonly observed among elderly subjects. Of note, the effects of ageing on the heart are minimal with respect to those in the main arteries [56].
As it was done in [29], we use the simulation of ageing to highlight the importance of accounting for heart-circulation interactions in cardiovascular modeling. In order to simulate the ageing of the vessel, we modify some relevant arterial parameters. We follow here the scheme proposed by [57], carrying out the appropriate manipulations to adapt this strategy to our model (single one-dimensional vessel). The parameter β, which is a surrogate of the arterial wall stiffness, the total arterial resistance and the total compliance are modified to represent "older" vessels. To take into account the contribution of the 1D model to the total resistance, the average pressure observed at the beginning of the aorta (P inlet ) is used as a marker for the impedance of the vessel. A target value (P target ) for this inlet pressure is calculated in order to induce the increasing pressure that is observed in Fig. 11 The important role of the cardiac model as an inlet condition for the distributed arterial model is highlighted comparing the output of the uncoupled aortic model (top left and right) and the fully coupled model (bottom left to right) settings, with a focus on the effects of ageing on the elastic arteries. On the left we depict the blood flow at the first element of the aortic model, in the center we show the aortic pressure (P ar ) at the same space element and at the bottom right we depict the pressure-volume loops in the left ventricle, for normal physiological conditions (black line) and for an arterial ageing conditions (blue line) normal ageing, and it is chosen as a 10% increase from the baseline inlet pressure. The peripheral resistance-R per in our RCR model-is updated conforming to where k is the ageing process iterative index. The parameter β is then modified as β new = 2.5β. The total arterial compliance is decreased as C tot,new = C tot /2, with the total compliance C tot being defined as The simulation of ageing is performed with the uncoupled aortic model and the fully coupled model settings.
The obtained results are presented in Fig. 11. In both simulation settings, we can observe numerous changes with respect to the baseline case that are typical of the ageing process: higher systolic peak, lower diastolic pressure and a small increase in wave propagation speed (the dicrotic notch is anticipated). However, the pressure curve obtained with the coupled model is more physiological. In particular, it is possible to see a trend towards the merging of the pressure systolic peak and the dicrotic peak, as it happens when the ageing process is advanced [56], and the dicrotic notch fades away. Indeed, without a cardiac model the inflow has to be imposed and therefore it does not adapt to the arterial conditions. In the ageing case with the uncoupled aortic model it is possible to observe a double reflection that does not reflect the arterial response to ageing. Furthermore, when using a coupled model, the blood flow is impacted itself by the ageing process. Our simulation shows a decreased peak compared to the baseline case. Finally, the use of a coupled heartcirculation model allows to compute the variation of the left ventricular pressure-volume loop as a result of the ageing process. This enables to monitor the quantitative effects of the arterial ageing on the cardiac function.

Conclusions
In this work, we present an original method to couple reduced-order blood flow circulation models and heart models through the point of view of energy balance both at the continuous and the discrete level. This is a difficult path as model reduction is often associated with the loss of energy balances for limit models. In fact, our coupled formulation is here proven to satisfy a full energy balance at the continuous level, with an additional consistent and controlled numerical dissipation on the fully discrete scheme. This allows to control the energy sources avoiding instabilities such as for instance uncontrolled back flows or perpetual motion cardiac engines. Essentially, our controlled coupling improves the modeling from the cardiac side and from the cardiovascular side. On the one hand, from the cardiac point of view we obtain better physiological signals with respect to the completely lumped heart-plus-Windkessel model proposed in [32]. On the other hand, we introduce a more physiological heart engine than a phenomenological time-varying elastance model acting as a flow generator in [29]. We believe that this paves the way for suitably investigating physiological aspects of heart-circulation coupling, as illustrated here with the dicrotic notch and the ageing process.