## Abstract

Complex congenital heart disease characterized by the underdevelopment of one ventricular chamber (single ventricle (SV) circulation) is normally treated with a three-stage surgical repair. This study aims at developing a multiscale computational framework able to couple a patient-specific three-dimensional finite-element model of the SV to a patient-specific lumped parameter (LP) model of the whole circulation, in a closed-loop fashion. A sequential approach was carried out: (i) cardiocirculatory parameters were estimated by using a fully LP model; (ii) ventricular material parameters and unloaded geometry were identified by means of the stand-alone, three-dimensional model of the SV; and (iii) the three-dimensional model of SV was coupled to the LP model of the circulation, thus closing the loop and creating a multiscale model. Once the patient-specific multiscale model was set using pre-operative clinical data, the virtual surgery was performed, and the post-operative conditions were simulated. This approach allows the analysis of local information on ventricular function as well as global parameters of the cardiovascular system. This methodology is generally applicable to patients suffering from SV disease for surgical planning at different stages of treatment. As an example, a clinical case from stage 1 to stage 2 is considered here.

## 1. Introduction

The hypoplastic left (or right) heart syndrome (HLHS or HRHS, respectively) is a complex congenital heart disease encompassing a wide range of cardiac anatomical abnormalities and is characterized by the underdevelopment of one ventricular chamber, which is invariably fatal during the first weeks of life without intervention. Normally, a three-stage surgical repair is performed to obtain separated pulmonary and systemic circulations (figure 1) [1,2]. The technique involves subsequent changes in the circulatory layout, up to a final configuration where the venous return passively flows into the lungs without the need for a pumping chamber, and the functional single ventricle (SV) provides the systemic blood flow. The SV heart is subjected to abrupt changes in the working conditions after each surgical stage as well as to an unusually high work that over time might cause ventricular dysfunction.

Over the past three decades, dramatic improvements in the treatment and outcomes for SV diseases have been achieved. Despite major improvements in terms of life expectations, current surgical strategies still remain palliative [3]. Great efforts should be paid to improve the surgical outcomes and consequently the patient quality of life and life expectation.

In this view, patient-specific computational models, able to account for different anatomies and for the associated secondary pathologies, can play an important role to help predict the surgical treatment outcomes. In addition, the improvements in the imaging techniques presently make it easier to build patient-specific models of cardiovascular haemodynamics for both healthy and pathological patients, e.g. the treatment of congenital heart diseases.

In recent years, several models have been proposed in the literature aimed at the evaluation of the surgical planning for congenital heart diseases [4–9]. These models, which describe the vascular area subjected to a surgical treatment in detail (three-dimensional geometries), suffer from a severe limitation owing to the inability to account for the interaction with the rest of the circulation. To overcome this important drawback, multiscale models have been proposed that couple the three-dimensional representation of the surgical site to a lumped parameter (LP) model [10,11] or a one-dimensional model of the circulatory system [12]. An LP model is a dynamic representation of the physics that neglects the spatial variation of parameters and variables, which are assumed to be uniform in each spatial compartment (i.e. zero-dimensional description). Conversely, one-dimensional models include the change of variables along the vessels and account for wave transmission effects within the vasculature. Therefore, an LP model gives rise to a set of ordinary differential equations describing the dynamics in time of the variables in each compartment, whereas one-dimensional models are described by a set of partial differential equations (axial coordinate and time). In the context of cardiovascular system modelling, the LP description is commonly applied to the major components of the system (e.g. cardiac chambers, heart valves and vascular compartments) to evaluate the global distributions of pressure, flow and blood volume.

Patient-specific multiscale models in a closed-loop fashion have been developed to simulate individual haemodynamics changes across those surgeries that significantly modify the circulatory layout [13–16].

Focusing on the detailed description of the SV behaviour, a few computational fluid-dynamic simulations of HLHS patient-specific cases have been recently performed to investigate the efficiency of the filling [17] and to evaluate the reliability of kinetic and viscous energy as clinical markers [18]. The diastolic function in HLHS patient models has been further investigated using personalized ventricular models by means of a fluid–structure analysis [19]. However, clinically measured ventricular inflow and outflow rates were used to prescribe the corresponding velocity boundary conditions in all these works. Hence, such three-dimensional models of SV cannot be applied in surgical planning to investigate the effect of different treatments, because different, unknown boundary conditions have to be prescribed in the post-operative conditions.

A closed-loop approach could solve this issue, because the coupled LP model of the circulation may provide the required new boundary conditions for the SV model.

Modelling SV defects and the change in working conditions owing to the surgical treatment is a challenging task because it requires accounting for different phenomena occurring at different scales. The behaviour of the SV at the organ level is strongly affected by the mechanisms occurring at the tissue level (sarcomere contraction). Moreover, the SV function results from the interaction with the whole circulatory system. Both the above-mentioned interactions call for a multiscale approach. Nevertheless, to the best of our knowledge, a multiscale model of SV cardiovascular system accounting for the three-dimensional anatomy and structure of the ventricle has not been proposed yet.

In this context, the aim of this study is to apply the multiscale computational framework developed by Krishnamurthy *et al.* [20] to couple, in a closed-loop fashion, a patient-specific three-dimensional finite-element (FE) model of the SV to patient-specific LP models of the pre- and post-operative circulations.

In addition to haemodynamic information (flows and pressures) through the circulatory system, the modelling of the three-dimensional structure of the SV allows the evaluation of local parameters as the myocardial stress and strain distributions, which are fundamental in the evaluation of the performance of the ventricle. Furthermore, global parameters such as pressure–volume (PV) loops and cardiac performance (i.e. cardiac output (CO) and ejection fraction) can also be calculated [21,22].

A key issue in the development of the model is the estimation of both the patient-specific SV properties and the circulatory parameters. To this aim, a sequential approach was carried out in this study. First, the circulatory parameters were estimated by using a closed-loop LP description of both SV and circulation (fully LP model). Then, the SV material parameters and the unloaded reference configuration were identified using the stand-alone three-dimensional FE model of SV. Finally, the fine tuning of these parameters was performed after coupling the three-dimensional model to the LP model of circulation, thus closing the loop. The circulatory parameter identification process was performed using the data clinically acquired before the surgery. Once pre-operative conditions were satisfactorily reproduced, the virtual surgery was performed, and the acute changes in the haemodynamics and the ventricle performance consequent to the surgery were predicted.

## 2. Material and methods

The model presented here relies on a multiscale strategy ranging from the three-dimensional FE model of the SV to the LP model of the circulatory system. The multiscale model comprises a significant number of parameters to be estimated for each patient-specific case, regarding both the material parameters (passive and active) of the myocardium and the resistances, the compliances and the inertances of the LP model. To this aim, and to reduce the computational demand, the sequential approach schematically illustrated in figure 2 was adopted consisting of three main steps. The parameter estimation was carried out by integrating the available clinical data, which were collected before surgery, and the literature information. In this study, the methodology was applied to a clinical case with a hypoplastic right ventricle.

### 2.1. Clinical pre-operative data

Clinical data for this study were collected at the Medical University of South Carolina (Charleston, SC) after institutional review board approval. Informed consent was obtained from the subject's parents. At the time of data collection, the patient was a three-month-old infant weighing 4.5 kg with a body surface area of 0.26 m^{2}. The baby was born with tricuspid atresia and pulmonary atresia with a severely hypoplastic right heart syndrome, unsuitable for a two ventricle circulation. The functional ventricle was the left one. On the fifth day of life, the patient underwent surgical placement of a right modified Blalock–Taussig shunt with a 4 mm PTFE (polytetrafluoroethylene) graft. His left pulmonary artery was severely stenotic and was augmented with a patch of bovine pericardium.

Clinical data were collected in terms of catheterization-derived pressure tracings, magnetic resonance (MR) flow tracings, MR-derived reconstruction of the SV and echocardiographic Doppler velocity tracings. Clinical measurements were performed prior to the stage 2 procedure. A pre-operative echocardiogram was performed under sedation and followed routine clinical protocols. MR was performed under general anaesthesia, and flow and volume measurements were obtained in different locations. Cardiac catheterization for pre-operative testing was performed under conscious sedation, using an ECG-gated, free-breathing cine-phase contrast velocity-encoded pulse sequence and commercially available cardiac analysis software. Haemodynamic assessment with intracardiac and vascular pressure measurements was obtained using a fluid-filled catheter system.

### 2.2. Development of the multiscale closed-loop patient-specific model

The patient-specific model was set up following three main steps.

In the first step (step 1, figure 2), the parameter estimation of the cardio-circulatory model was performed by means of a fully LP closed-loop model (figure 3, described in detail in the §2.3.1). The patient-specific parameters of the cardio-circulatory model were obtained through an optimization process described elsewhere [23]. Briefly, a Bayesian estimation of the parameters was performed. Adaptive Markov chain Monte Carlo was employed to obtain the distributions of the model parameters. Then, Nelder–Mead hill-climb optimization was performed from the parameter set that was found to maximize the posterior distribution during the Markov chain Monte Carlo iterations. Optimization was performed to find the patient-specific parameters that best replicated the clinical pre-operatory haemodynamics of the patient in terms of mean, maximum and minimum values of pressures, flows and cardiac volumes. In this phase, the use of a fully LP model to estimate patient-specific circulatory parameters was mandatory to dramatically reduce the computational effort.

In the second step (step 2, figure 2), the three-dimensional ventricular geometry of the patient was reconstructed from MR images, and the myocardial parameters (passive and active models) were estimated considering the stand-alone three-dimensional FE model of the SV (uncoupled from the circulatory model). First, the unloaded reference state and the myocardial passive properties were identified by means of the iterative method previously developed by Krishnamurthy *et al.* [20]. Then, the three-dimensional FE model was inflated to three different volumes and isovolumic contractions were simulated, iteratively tuning the active myocardial parameters. Both tuning processes made use, as targets, of the active and passive ventricular elastances obtained from the fully LP optimized model. The tuning strategies are described in §2.3.

Finally (step 3, figure 2), proper haemodynamic boundary conditions were provided to the three-dimensional FE model of SV through its coupling to the LP model of univentricular circulation; a multiscale closed-loop cardio-circulatory model was then created. The results previously obtained from the fully LP model were used as the initial conditions for the multiscale model. At this point, a further fine-tuning of active parameters was carried out to ensure that the multiscale model well replicated the clinical data of the patient.

In both steps 2 and 3, the three-dimensional model was subjected to proper kinematic boundary conditions to prevent displacement and rotation of the base along the long axis.

The optimization process and the fully LP closed-loop circulatory model were carried out in Matlab (The MathWorks, Inc., Natick, MA), whereas the FE simulations were performed using Continuity 6.4 (National Biomedical Computation Resource, University of California, San Diego, CA).

### 2.3. Components of the multiscale closed-loop patient-specific model

In the following sections, a detailed description of the various components of the LP circulatory model and the three-dimensional FE SV model, composing the multiscale model, is provided.

#### 2.3.1. Circulatory model

Two LP models of the circulation were developed corresponding to stage 1 and stage 2 of Fontan surgical procedures. The cardio-circulatory layout (figure 3) comprised four main blocks to model the single heart, the systemic upper (UB) and lower (LB) body circulations and the lungs. For the sake of simplicity and in the absence of any clinical information, coronary circulation, gravity and respiratory effects were neglected in the present case study.

The single atrium was described as a time-varying elastance which may actively contract according to a sinusoidal activation function. Both the atrioventricular and outflow valves were described by nonlinear diodes to enforce unidirectional flow. Valve inertial effects were directly implemented with an inertance for the atrioventricular valve and included in the downstream arterial parameters for the outflow valve. Valve parameters were calculated as a function of the measured valve area [24].

Each peripheral vascular region was modelled by means of an RLC–R–CR block describing arteries, capillary bed and veins. For stage 1, a resistive block mimicking the systemic-to-pulmonary shunt was added, including linear and nonlinear components accounting for both distributed and localized energy dissipations. Shunt parameters were expressed as functions of the shunt diameter as described in Migliavacca *et al*. [25]. The shunt model was removed when stage 2 was simulated (see §2.4).

To perform the parameter identification (previously described as step 1 of figure 2), the stage 1 LP circulatory model was coupled with an LP model of the SV (figure 3, red box (*a*)). This simplified LP model of the SV consisted of a time-varying nonlinear elastance (accounting for the active and passive properties of the ventricle) in series with a resistance mimicking the viscous behaviour of the myocardium [13,16]. The myocardium contraction was ruled by a sinusoidal activation function. The use of an LP of the SV model was mandatory in this step to obtain parameter estimation with a reasonable computational effort. After the circulatory parameters were correctly identified, the LP circulatory network (i.e. atrial and vascular models) was coupled to the three-dimensional FE model of the SV to create the multiscale model (figure 3, red box (*b*)). Lastly, the LP circulatory layout was changed as illustrated in figure 3 (orange colour) to simulate the post-operative state. For a detailed description of the general equations governing the LP circulatory model, the reader is referred to appendix A.

#### 2.3.2. Three-dimensional finite-element model of the single ventricle

A model of the left ventricular anatomy was reconstructed from MR images of the investigated paediatric patient, suffering by hypoplastic right heart syndrome. The manual segmentation was performed based on 50 slices equally spaced of 1 mm with respect to the ventricular long axis and with a spatial in-plane resolution of 0.25 × 0.25 mm. Because the *in vivo* MR images were acquired at the end-diastole (ED), the reconstructed geometry was further processed (see §§2.3.3). Indeed, to correctly compute stresses, an unloaded reference state must be identified, defined as the state at which the ventricle cavity pressure, in a passive state, is zero.

The domain discretization of the ED geometry was performed using hexahedral cubic-Hermite FE mesh. Geometries represented by means of Hermite interpolation guarantee C^{1} continuity on the boundary between adjacent elements, as the degrees of freedom of this mesh are both nodal coordinate values and nodal derivatives. This description yields a smooth description of the shape, which represents a great advantage when dealing with soft tissue simulations.

First, the endocardial and epicardial surfaces of the SV at ED were identified and manually segmented (figure 4*a*). Second, a simplified two-dimensional cubic-Hermite template mesh (80 nodes and 64 elements) was defined in a prolate spheroidal coordinates system consisting of two nested surfaces representing the endocardial and epicardial surfaces. The template was then fitted in the radial coordinate to the reconstructed data (figure 4*b,c*). For a detailed description of the method the reader is referred to reference [26]. The obtained surface meshes were then converted into three-dimensional hexahedral cubic-Hermite elements in rectangular Cartesian coordinates. In order to better define the transmural fibre orientation gradient, the mesh was refined across the myocardial wall thickness to a final resolution of 120 nodes and 64 elements (figure 4*d*).

The three-dimensional FE model also included a myofibre and sheet architecture. Fibres were assumed to be tangential to the endocardial and epicardial surfaces, so that the material anisotropy of the ventricular wall was defined by referring to an orthogonal system of local material coordinates, having one axis aligned with the muscle fibre direction. The fibre angle at any point in the model was given by an interpolation of the fibre field parameters defined at the same nodal positions used to define the ventricular geometry. The measured transmural distribution of fibre and sheet angles was approximated using trilinear Lagrange basis functions.

As specific information on the patient-specific fibre orientation was not available, fibre architecture was assumed as in physiological adult ventricles. Because the functional ventricle of the considered patient was the left ventricle, a fibre orientation ranging from −60° at the endocardium to +60° at the epicardium (with respect to the circumferential direction) was assumed [27].

#### 2.3.3. Passive properties and unloaded geometry

The passive material properties of the myocardium were described using the transversely isotropic form of the Holzapfel constitutive model [28]. If the myocardium is described as an incompressible medium with only one family of fibres perfectly aligned along their mean direction, the strain energy function includes two exponential terms, accounting for the isotropic matrix and the anisotropy of the fibres. For a detailed description of the constitutive equation the reader is referred to appendix B. Two material parameters had to be estimated to reproduce a patient-specific behaviour of passive ventricle.

The myocardial passive material parameters were estimated together with the identification of the unloaded geometry. To this aim, an iterative method previously developed by Krishnamurthy *et al.* [20] was adopted. This unloading algorithm (step 2 in figure 2) was able to estimate the unloaded geometry on the basis of the ED ventricular model, if the ED pressure (measured) and the unloaded volume (assumed) were known, whereas the passive material parameters were iteratively tuned. As the initial guess, the material parameters of the model were chosen based on literature experimental biaxial tests performed on ex vivo canine myocardial tissue [29,30]. Starting from this information, the passive material parameters were manually tuned, until an unloaded geometry that matched the unloaded volume was obtained. This unloaded geometry, when inflated to ED pressure, led to a loaded geometry that well approximated the original ED fitted mesh.

The unloading algorithm was based on the multiplicative decomposition of the deformation gradient tensor, as described in tissue growth modelling [31]. The iterative estimation scheme aimed at finding the unloaded reference geometry that minimized the difference between the measured ED geometry and the computed geometry when the unloaded model was inflated to the measured ED pressure, given nonlinear myocardial material properties. In the first step, the fitted ED model was chosen to be the first guess unloaded state and inflated to the measured ED pressure in the patient. The inverse of the resulting deformation tensor (computed between the inflated geometry and the fitted ED geometry) was then applied to the measured ED geometry to compute a new unloaded geometry. If the associated volume deviated from the assumed unloaded volume (*V*_{SV0}), the material stiffness was changed (e.g. if the FE model unloaded volume was smaller than the assumed volume, the material stiffness was increased). This sequence was then iterated until a suitable unloaded geometry was obtained with a cavity volume equal to *V*_{SV0}.

In the work of Krishnamurthy *et al.* [20], the estimation of the unloaded volume was based on an empirical relation [32], as a function of the ED pressure and volume, which was found to fit the measurements in human adults and animal hearts in healthy and disease conditions. However, this relation could not be applied to paediatric patients because it provided unrealistic results (unloaded volumes always larger than end-systolic volume). Hence, in this work, the unloaded volume resulting from the circulatory layout parameters optimization (step 1) was adopted as *V*_{SV0}.

Once the unloaded geometry was obtained and the appropriated passive material properties were identified, the ED geometry measured from the MR images was obtained by inflating the unloaded geometry up to the ED pressure.

#### 2.3.4. Active material model

The model proposed by Lumens *et al.* [33] was adopted to describe the active contractility properties of the myocardium. This phenomenological representation of sarcomere contraction simulates experiments on isolated rat cardiac muscle. For a detailed description of the active model, the reader is referred to appendix B.

Briefly, the sarcomere was modelled as a passive element in parallel with a series combination of a contractile element and an elastic element. The time-dependent law of the contractile element depends on two state variables, i.e. the length of the contractile element and the mechanical activation. The mechanical activation function represents the complex phenomena occurring at the intracellular level, whereas the sarcomere length is obtained from the myofibre strain. The length of the series elastic element and the myofibre force–velocity relation linearly determines the time derivative of the sarcomere length. The time derivative of the mechanical activation separately accounts for the rise and decay of the mechanical activation. The parameters of the myocardium active model to be estimated are the active stress scaling coefficient and the scaling factors of the activation time-function. A suitable patient-specific value for the first parameter was assessed by manual tuning, and use of the active ventricular elastance obtained from the fully LP optimized model (see step 2 in figure 2). Namely, the stand-alone three-dimensional FE model of the SV was inflated up to three different volumes. Isovolumic contractions were simulated, and the ventricular pressure was evaluated at maximum contraction. The active stress scaling coefficient was iteratively tuned until the obtained pressures matched those deduced from the fully LP model. Finally, starting from the activation function obtained from the fully LP optimized model, the scaling factors of the contraction time-function were finely calibrated until the response of the multiscale closed-loop model well reproduced the clinical data.

### 2.4. Simulation of the surgical procedure

Once the pre-operatory conditions of the clinical case were satisfactory simulated with the multiscale model, a virtual surgery was performed by changing the layout of the LP circulatory model. Figure 3 shows the difference between stage 1 and stage 2 circulatory layouts: the shunt was removed, the UB circulation was disconnected from the atrium and connected to the lungs. In the absence of specific information, the parameters of the circulatory model were assumed to be unchanged from the pre-operative situation, thus simulating an acute change. Similarly, the myocardial properties of the three-dimensional SV model remained unchanged.

Simulations results were represented in the form of global circulatory outputs (both mean values and time tracings) and local effects in the cardiac wall (myofibre stress and strain at different time points of the cardiac cycle).

The predicted post-operative conditions were compared with the pre-operative ones to assess haemodynamic effects of the performed surgery.

## 3. Results and discussion

### 3.1. Set-up of the patient-specific pre-operative model

The results of the simulation of the patient pre-operatory conditions after the optimization process are summarized in table 1. Results from the fully LP model and the multiscale model were reported against the pre-operative clinical data. An overall good agreement was obtained for both models in describing the main clinical quantities, with an error lower than 5%. A proper CO indicates a satisfactory description of ventricular–arterial and venous–atrial couplings. Moreover, the correct repartition of flows to the systemic (*Q*_{ub} and *Q*_{lb}) and pulmonary circulations (*Q*_{p}) ensures a good modelling of the systemic-to-pulmonary shunt.

Focusing on the multiscale model results, figure 5 shows the pressure–volume loop (*a*) and pressure tracings (*b*) of the single ventricle (*P*_{sv}) and of the aorta (*P*_{ao}) for three consecutive cardiac cycles. It is worth observing that a steady behaviour was obtained after only two cycles (difference lower than 1% between subsequent cycles). The use of the initial values provided by the fully LP model allowed us to limit the simulation to a few cardiac cycles, with an effective cost reduction in computation time.

It is also worth noting that although the active material parameters of the stand-alone three-dimensional SV model were adjusted on the basis of the elastances obtained after the fully LP model optimization, a finer tuning was still required to match the patient behaviour when the multiscale model was used. This was mainly due to the different description of the myocardium contractility, in particular to the effects of the velocity of contraction. Indeed, while the LP model of SV in the fully LP model accounted for a constant viscosity of the myocardium acting when the SV volume changes, the viscous term included in the active model of the three-dimensional model relied on a more complex description based on sarcomere mechanics [33]. However, only a few iterations were needed to obtain a good prediction of the patient pre-operatory behaviour, thus limiting the overall computational cost.

### 3.2. Simulation of stage 1 to stage 2 surgery

The changes in the circulatory layout, simulating the surgery between stage 1 and stage 2, led to an increase in the afterload and to a decrease of the preload of the SV. Figure 6 shows the pressure tracings (*P*_{sv}, *P*_{ao}) of six consecutive cardiac cycles after the circulatory model was changed to simulate the surgery. Simulation of six cardiac cycles was necessary to reach a steady state (difference in pressure lower than 1% between two subsequent cycles).

Figure 7 shows the pre-operatory and post-operatory pressure–volume loops. As a consequence of the surgery, it is interesting to observe a significant decrease of the stroke volume (from 13.7 to 10.4 ml) with a slight increase in end-systolic volume together with larger decrease of ED volume.

The comparison of the multiscale haemodynamic results for stage 1 and stage 2 conditions is reported in table 1 (last two columns). It is worth noting that the post-operative model predictions were not validated, because post-operative clinical data were not available owing to the patient's good and stable condition. A decrease in the CO of about 25%, together with an increase of the ventricular peak pressure (17%), can be noted. Mean aortic pressure increases (from 53.9 to 74.1 mmHg); similar increments also occurred in both systolic and diastolic values. Conversely, the atrial pressure decreased of about 40%. These significant changes were expected with ventricle volume unloading and afterload increasing. In stage 1 circulation, the total resistance was composed of three resistances in parallel: UB and LB systemic vascular resistances and the series of shunt and pulmonary vascular resistances (PVR). After stage 2 surgery, the shunt resistance was removed, whereas the UB is connected in series with the PVR. As a consequence, the total afterload of the heart increased post-operatively. The ventricular unloading represented a beneficial change, because the ventricle operated at a lower ED volume, thus reducing the risk of critical conditions that could lead to sarcomere decompensation.

The stress and strain distributions in the myocardial wall were evaluated for both the stage 1 and stage 2 conditions. Figure 8 shows the stress pattern for both the pre-operative and post-operative conditions at ED (upper panels) and at systolic peak (lower panels). Concerning ED, a reduction of the maximum stress in the myocardium was observed from stage 1 (5 kPa) to stage 2 (3.9 kPa) as well as in the overall values. This reduction of myocardial stress was a consequence of the change of the layout of the circulation, which leads to the unloading of the ventricle. As expected, the stress pattern distribution was preserved in the two cases. Indeed, the simulation of the stage 2 surgical procedure represented an acute condition, thus the remodelling process associated with chronic condition was neglected. With respect to the systolic peak, an opposite trend was observed. Indeed, a higher systolic peak pressure was found in stage 2 with respect to stage 1 as shown in figure 7.

The end-diastolic strains in the fibre direction are depicted in figure 9, at the stage 1 and stage 2 end-diastole. An overall reduction of the fibre strains in the stage 2 was found with respect to stage 1. In stage 1, significant portions of the myocardial wall were subjected to fibre strains ranging from 20% to 25%. These values are known to be rather critical because they are very close to sarcomere decompensation. Both for stage 1 and stage 2 conditions, the time profile of the myocardial strains computed at different locations (points at epicardial and endocardium middle cross section and at the apex) are shown in figure 9. A reduction of maximum strain in stage 2 was observed in all locations. A lower strain variation from ED to end-systole was observed in stage 2 (endocardium: 16%; epicardium: 14.5%; apex: 20.7%) with respect to stage 1 (endocardium: 22%; epicardium: 19%; apex: 23%).

### 3.3. Possible use of the model for clinical support

The multiscale model presented might be useful in the evaluation of the acute post-operative outcomes. Indeed, the surgical procedure leads to an abrupt alteration of the working conditions of the SV which are typically hardly predictable. In addition, the timing of the surgical treatment is critical and based on the clinical experience of the surgeons. In this view, the multiscale model can be used to simulate a number of different surgical scenarios, thus providing important information to assist the surgeon in the decision-making process. Moreover, this information can help customize the post-operative pharmacological treatment.

The systematic application of this methodology to a number of SV patients could provide further insights into material properties and fibre orientation data for this pathology, still lacking in the literature. Indeed, the material parameters and fibre orientation could be calibrated, so that the kinematic outcomes of the multiscale model (i.e. myocardial wall displacement and deformation) reproduce the four-dimensional displacement field of the patient. In addition, the study of the two different groups of SV patients (HLHS and HRHS) could highlight both inter- and intragroup variability, thus improving the understanding of such a complex pathology.

## 4. Limitations

The proposed computational framework represents a robust tool to couple a patient-specific three-dimensional model of SV and an LP model of the circulatory network. The framework easily enables the simulation of different circulatory models to describe several surgical treatment scenarios. Nevertheless, the model has a number of limitations.

First, the main limitation of this study is that the model results were not validated against clinical data. The model parameters were identified to replicate the pre-operative patient haemodynamic quantities. Concerning the pre-operative state, a validation process will be performed as soon as more detailed clinical data (e.g. four-dimensional displacement field of the myocardium) are available for stage 1 patients. Regarding the post-operative results, no data were collected for the investigated patient, as they would require invasive and unnecessary clinical exams. In addition, the simulations of the surgical treatment were performed assuming unchanged material and circulatory model parameters. However, no quantitative information is presently available to describe possible active changes or adaptations of the cardiovascular system in this kind of paediatric patients.

Second, owing to the lack of patient-specific measurement, a fibre architecture consistent with the literature physiological data was assumed. Fibre orientation in pathological conditions may be significantly different with respect to the physiological case. Moreover, the estimation of the material anisotropy may be incorrect. To this issue, the combination of the three-dimensional cine MR images measurement with the already available clinical data can provide important information. Indeed, the four-dimensional displacement field of the myocardium might be used to validate our assumptions.

## 5. Conclusion

In this paper, a multiscale approach used to simulate patient-specific biomechanics in the case of SV defect was presented. The modelling process integrated clinical data, including MR imaging, catheterization and echocardiography, to develop a patient-specific model able to account for the main features of the patient physiology. Our approach involved distinct steps: optimization of patient-specific parameters of both the three-dimensional SV and the LP circulatory model, creation of the multiscale pre-operative model and modification of the circulatory layout for post-operative simulations.

The methodology proposed was then applied to model a single patient-specific stage 2 surgery. Nevertheless, the methodology is generally applicable to further patient-specific cases suffering from both hypoplastic left heart syndrome and hypoplastic right heart syndrome. Virtual surgical planning at different stages of the treatment can also be easily simulated by changing the layout of the circulatory model.

## Funding statement

This study was supported by a grant from the Fondation Leducq, Paris, France. The authors gratefully acknowledge use of the software Continuity 6.4 granted by the National Biomedical Computation Resource (NIH grant no. P41 GM103426).

## Acknowledgements

The authors acknowledge the use of the optimization code granted by Dr Daniele Schiavazzi from University of California, San Diego, CA, USA and Dr Alessia Baretta from Politecnico di Milano. The authors thank Dr Chiara Corsini from Politecnico di Milano for the useful discussions during the preparation of the manuscript.

## Appendix A. Equations of the circulatory models

The general equations used in the LP description of the circulation are presented here.

All the circulatory blocks (UB and LB circulations, right and left pulmonary circulations) were described by means of RLC–R–CR blocks. The following equations (one algebraic and two differential) were used to describe the blood pressure and volume.
A 1
A 2
A 3
In the equations, *R*, *C* and *L* refer, respectively, to the generic resistance, compliance and inertance of each vascular block. The terms *P*_{up}, *P*_{dw}, *Q*_{up} and *Q*_{dw} represent the generic pressures (*P*) and volumes (*Q*) upstream and downstream from each component.

The aortic valve was described by means of diode with nonlinear resistance *K*_{ao}. The pressure drop across the valve was calculated based on
A 4
where *P*_{sv} is the SV pressure, *P*_{ao} and *Q*_{ao} represent the pressure downstream and the flow across the valve. *R*_{myo} is a term accounting for the viscous behaviour of the myocardium and *K*_{ao} in the nonlinear resistance of the valve.

The flow across the valve was then ruled by the following equations A 5

The atrioventricular valve was modelled as a diode with nonlinear resistance in series with an inertance. The pressure drop across the valve is
A 6
where *P*_{sa} is the pressure in the single atrium, *L*_{av} and *K*_{av} are the atrioventricular inertance and nonlinear resistance and *Q*_{av} is the atrioventricular flow.

The flow across the valve was defined as follows A 7 Single-atrium and single-ventricle were described by means of time-varying elastances accounting for the active and the passive behaviour of the myocardium. The pressure generated by a generic chamber was defined as A 8

where ch = sv, sa indicates the single-ventricle and the single-atrium. The viscous term *R*_{myo} was considered only for the ventricle. Ventricular active pressure–volume relationship was considered nonlinear, whereas a linear relation was considered for the atrium. Both atrial and ventricular passive pressure–volume relationships were described with a nonlinear relation.

The following equations were used
A 9
A 10
A 11
A 12
where *V*_{sv0} and *V*_{sa0} are the unstressed volumes of each cardiac chamber (i.e. the volume at zero pressure), *a*, *b*, *c*_{sv} and *d*_{sv} are model parameters specific to the single-ventricle, whereas *C*_{sa}, *c*_{sa} and *d*_{sa} are model parameters specific to the single atrium.

The term *A*_{ch}(*t*) in equation A 8 ranges from 0 during diastole and 1 at the end of systole. It describes the excitation–relaxation pattern of the sarcomeres and is defined as
A 13
A 14
where *t*_{s} and *t*_{sa} are atrial and ventricular systole durations.

Model parameters for the patient-specific case are available upon request to the authors.

## Appendix B. Models of passive and active mechanics

The equations modelling the passive and active mechanics are described here.

The transversely isotropic form of the constitutive model developed by Holzapfel & Ogden [28] was assumed to model the passive properties of the resting myocardium. In this model, the anisotropy in the fibre and cross-fibre directions of the myocardium has been accounted using separate exponential terms with different exponents. Equation (B 1) gives the equation for the strain energy
B 1
where *a*, *b*, *a*_{f}, *b*_{f} are material constants, the *a* parameters having dimension of stress, whereas the *b* parameters are dimensionless. *I*_{1} is the first invariant of the right Cauchy–Green strain tensor (isotropic term), *I*_{4f} are the components of the right Cauchy–Green strain tensor in the fibre direction (transversely isotropic term).

The active model from Lumens *et al*. [33] was used to model the active behaviour of the myocardial fibres. In this model, the sarcomere has been modelled as a passive element in parallel with a series of a contractile element and an elastic element. The active stress (*σ*_{f,act}) added to the Piola–Kirchoff stress tensor to account for the contribution of fibre contraction, was defined as
B 2
where *s*_{act} is a factor scaling the active myofibre stress, *C* the mechanical activation, *L*_{sc} the contractile element length, *L*_{sc0} the contractile element length with zero load, *L*_{sc} the sarcomere length and *L*_{sc,iso} the length of isometrically stressed elastic element. The active stress depends on two state variables (i.e. the sarcomere length *L*_{sc} and the mechanical activation *C*). The time derivative of *L*_{sc} can be written as
B 3
The time derivative of *C* is
B 4
where *τ*_{R} and *τ*_{D} are scaling rise and decay time, respectively, *t* and *C*_{rest} represent time and diastolic resting level of activation, respectively, and *C*_{L}, *F*_{R} and *T* are functions describing the increase of activation with the sarcomere length, the rise of mechanical activation and the decrease of activation duration with decrease of sarcomere length, respectively.

## Footnotes

↵† Modeling of Congenital Hearts Alliance (MOCHA) Group: Andrew Taylor, Alessandro Giardini, Sachin Khambadkone, Silvia Schievano, Marc de Leval and T. -Y. Hsia (Institute of Cardiovascular Science, UCL, London, UK); Edward Bove and Adam Dorfman (University of Michigan, Ann Arbor, MI, USA); G. Hamilton Baker and Anthony Hlavacek (Medical University of South Carolina, Charleston, SC, USA); Francesco Migliavacca, Giancarlo Pennati and Gabriele Dubini (Politecnico di Milano, Milan, Italy); Alison Marsden (University of California, San Diego, CA, USA); Jeffrey Feinstein (Stanford University, Stanford, CA, USA); Irene Vignon-Clementel (INRIA, Paris, France); Richard Figliola and John McGregor (Clemson University, Clemson, SC, USA).

One contribution of 11 to a theme issue ‘Multiscale modelling in biomechanics: theoretical, computational and translational challenges’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.