Quantitative Study of the Effect of Tissue Microstructure on Contraction in a Computational Model of Rat Left Ventricle

Tissue microstructure, in particular the alignment of myocytes (fibre direction) and their lateral organisation into sheets, is fundamental to cardiac function. We studied the effect of microstructure on contraction in a computational model of rat left ventricular electromechanics. Different fibre models, globally rule-based or locally optimised to DT-MRI data, were compared, in order to understand whether a subject-specific fibre model would enhance the predictive power of our model with respect to the global ones. We also studied the impact of sheets on ventricular deformation by comparing: (a) a transversely isotropic versus an orthotropic material law and (b) a linear model with a bimodal model of sheet transmural variation. We estimated ejection fraction, wall thickening and base-to-apex shortening and compared them with measures from cine-MRI. We also evaluated Lagrangian strains as local metrics of cardiac deformation. Our results show that the subject-specific fibre model provides little improvement in the metric predictions with respect to global fibre models while material orthotropy allows closer agreement with measures than transverse isotropy. Nonetheless, the impact of sheets in our model is smaller than that of fibres. We conclude that further investigation of the modelling of sheet dynamics is necessary to fully understand the impact of tissue structure on cardiac deformation.


Introduction
Cardiac tissue is arranged into a complex architecture. Due to their elongated shape, myocytes are highly anisotropic and their local orientation in cardiac tissue is usually called the fibre direction [1]. Furthermore, the extracellular connective tissue matrix that surrounds myocytes groups them into layers 4-6 cells thick [2], separated by cleavage planes, forming a branching laminar architecture. Such layers are commonly referred to as sheets. Although there is general agreement about the local alignment of cells and their organisation into sheets, disagreement remains on their exact characteristics in beating hearts, their alteration due to pathology and in general their degree of contribution to cardiac function [3]. Tissue microstructure contributes to the heterogeneity of tissue conductivity and therefore directly affects electrical wave propagation, with dramatic consequences in the case of pathology. Regarding mechanics, the importance of taking into account fibre and sheet orientation in the definition of passive and active behaviour of cardiac tissue has already been highlighted [4]. It is along the fibre direction that electrical excitation propagates fastest and that the principal component of contraction takes place. When included in an electromechanical model, fibres are often mathematically prescribed [5]. Recent advances in imaging technology have made it possible to measure local microstructure from data such as ex vivo diffusion tensor MRI (DT-MRI). The primary eigenvector of the diffusion tensor has been shown to correlate with the fibre direction in the tissue [6]. Electromechanical models based on DT-MRI data suggest that fibres are a dynamical system in cardiac structure, undergoing adaptation that minimises shear strains, behaviour that can have important consequences in tissue remodelling in pathology [7]. The importance of local myocyte orientation in cardiac function, both in health and disease, has promoted considerable effort in recent years towards subjectspecific modelling of cardiac electromechanics [8], [9], [10]. Pathological conditions such as dyssynchrony [11] have been investigated, with realistic ventricular to whole-organ geometries, combined with fibre orientations obtained from imaging [10]. The results are promising, but it is an open question as to whether current electromechanical models are becoming a useful support tool in clinical decision-making.
The laminar structure of myocardial tissue is also important for cardiac electromechanics. It has been established that sheets undergo shortening in the fibre direction, extension in the insheet/cross-fibre plane, thickening and shear, whose regional variability and transmural distribution are important for the distribution of strain and stress throughout the ventricle during cardiac motion [12]. More recently, changes in sheet intersection angles have been identified to occur during the cycle of contraction and relaxation [13]. These sheet dynamics are thought to form one of the underlying mechanisms to explain wallthickening during contraction [14]. Sheets have been modelled as continuously varying transmurally [2], although more recent experimental findings suggest that the transmural distribution of sheets is more complex and discontinuous [13], see for example Figure 1(f). Sheet orientations directly obtained from histological data, usually from a small number of ventricular locations only, have been employed in computational studies of the effect of tissue microstructure on electrophysiology [15], but in general the importance of sheets has not been fully acknowledged in the cardiac modelling community in spite of the experimental evidence.
Two motivations underlie our study. First, we aim to understand whether the use of subject-specific fibres would enhance the predictive power of our model, or whether it would be as informative as simple globally prescribed fibres within our modelling framework. To this purpose, globally defined rule-based fibres were compared with a subject-specific model of fibres that is locally optimised in order to minimise the discrepancy with DT-MRI data [5].
Our second motivation is the exploration of the often underestimated importance of sheet dynamics in the common electromechanical models. We want to test the impact of sheet orientation in the mechanics of the left ventricle. This was achieved by (a) comparing transverse isotropic with orthotropic material laws for cardiac tissue and (b) under the orthotropic material properties, testing two global models of transmural variation of the sheet orientation, namely a linear and a bimodal model. The aim of the latter comparison was to see whether the latest findings regarding the discontinuous multiple-population distribution of sheets [13] (mimicked by the bimodal model) can improve the accuracy of the electromechanical model, as opposed to the standard smooth variation of sheet orientation [2] (the linear model).
In order to quantify the effect of the different models of fibre and sheet orientation on contraction, global and local metrics were considered. The global metrics are ejection fraction, wall thickening, and shortening in the base-to-apex direction. We compared their predicted values with those obtained from rat cine-MRI data. The local metrics considered are the Lagrangian circumferential strain, as a measure of contraction, circumferential-radial shear strain, as a measure of the relative twist of sub-endocardial and subepicardial layers, and circumferential-longitudinal shear strain, as a measure of torsion [16].
A preliminary study of the effect of fibre direction and material orthotropy on contraction had been performed with a simplified electromechanical model (without pressure boundary conditions) [17]. The current study represents a significant extension in our model development, to further the study of the impact of sheets and the analysis of the local effect of tissue structure on deformation. The contribution of the study is twofold. First, we show that the fibre model affects significantly the prediction of the metrics considered. Therefore, it is important to choose the most suitable fibre setting. On the other hand, our results question the efficacy of using personalised myocyte orientations when compared to the limitations of state-of-the-art electromechanical models, since they do not seem to provide substantial improvements, both in terms of global and local measures, compared to simpler global rule-based fibre definitions, at least for the healthy left ventricle of rat. Secondly, we highlight the need for further development of models of sheets. On the one hand, the introduction of sheet orientation in the model by switching from transverse isotropy to orthotropy shows a small but consistent improvement in the metrics predictions. On the other hand, we found that, although experimentally more realistic, the bimodal model for sheet orientation did not improve and in some cases counteracted the effect of introducing material orthotropy. We believe that such an important structural change should affect more evidently the predicted deformation patterns. Given the results presented in this study, we conclude that our current electromechanical modelling approach may not be able to successfully make use of subject-specific tissue structure, with all its additional complexity and computational cost, unless a better understanding of the biological roles, and the computational modelling of sheets and their mechanical properties is achieved.

Materials and Methods
No new animal experiment was required for the results presented in this paper. The images used come from a dataset previously published and for which the ethics statement was duly reported [18]. However, with respect to those experiments, the authors confirm that all animal work was conducted in accordance with Schedule 1 of the UK Home Office Guidance on the Operation of Animals (Scientific Procedures) Act of 1986, and was approved by the Oxford University ethical review board.

MRI dataset
In vivo cine-MRI and subsequent ex vivo DT-MRI of rat hearts were acquired on a 9.4T (400 MHz) MR system (Agilent, Santa Clara, CA, USA). The cine-MRI dataset used a 2D multi-frame gradient echo sequence, and consisted of a stack of eight contiguous short-axis slices covering the heart from apex to the atrioventricular separation. Each short-axis slice comprised of 32-36 frames acquired throughout the cardiac cycle with an in-plane resolution of 100x100mm and a slice thickness of 1:5 mm. A shortaxis view in end-diastole is shown in Figure 1(c). The DT-MRI images where subsequently acquired on excised rat hearts, fixed either in resting state or in contracture and embedded in agarose, at an isotropic resolution of 100mm as described in [18].

Geometrical mesh for rat left ventricle
We approximate the left ventricle with a thick-walled truncated ellipsoid. A first coarse tetrahedral mesh was generated with semiaxes and wall thickness estimated from the end-diastolic frame of one of the cine-MRI rat scans described (Figure 1(c)), which we took as the reference mesh (Figures 1(a) -1(b)). In order to estimate the shape of a mesh representative of the left ventricle in absence of blood pressure (called pre-P mesh), an iterative approach was applied. The reference mesh was shrunk by applying a uniform scaling matrix to give the pre-P mesh, which was then allowed to expand passively under a constant pressure of {1 KPa applied to the endocardial surface to give the post-P mesh (the negative sign means that the pressure is applied inwards along the normal to the surface). This step was repeated for several iterations adjusting the scaling factor until we obtained an optimal pre-P mesh producing a post-P mechanics mesh close to the reference mesh. Details about the pre-P and post-P meshes are given in Table 1. The mechanics mesh for the electromechanical simulations is the pre-P mesh resulting from the optimisation performed, and has quadratic elements. The mesh for the electrophysiological component is instead obtained by refining the post-P mechanics mesh and has approximately 1:6 : 10 4 linear tetrahedral elements with a mean edge length of 5:7 : 10 {2 cm. All meshes were generated with the open-source 3D finite element grid generator Gmsh [19].

Tissue structure definition
Tissue structure is defined at the centroid of each element of the mechanical mesh by assigning fibre, sheet and sheet-normal directions. This orthonormal local system of coordinates is assigned according to a well-established approach [2]. First, the orthonormal system of material coordinates (radial u, circumferential v and longitudinal w directions) for the centroid of a given element is identified. A first rotation about the radial direction u is performed, that transforms the vectors v and w into v a and w a . The rotation angle a is the helix angle that defines the fibre direction f~v a . The second step consists of a rotation about the fibre direction f of the angle b, the sheet angle. The vector u is rotated to u b that represents the sheet direction s, while the vector w a is rotated into w ab that gives the sheet-normal direction n. Different models for the fibre and sheet orientation are introduced in this tissue structure generation framework by modelling the transmural variation of the helix and sheet angle as described below.
Models of fibre direction. A set of rule-based globally defined fibres were generated by assigning to each element of the mesh a value of the helix angle a, defined as the function a~a(R,n,d) and calculated by eqn. (1) a~R Ã sgn(1{2d)j1{2dj n ð1Þ where sgn stands for the sign function, d is the normalised transmural depth, varying from 0 to 1 from endocardium to epicardium, R is the transmural rotation, i.e. the maximum absolute value of a allowed by the model, which will span between {R at epicardium and zR and endocardium. Finally, the exponent n modulates the transmural variation of a, so for n~1 it will be linear, n~3 cubic. This is a standard approach for globally prescribed fibres introduced by Streeter et al. [1]. Four values of R were tested, R~50, 60,70, 80 degrees, spanning the range of values found experimentally [5]. For each value of R, four values of the exponent n were considered, n~0:5,1, 2, 3, for a total of 16 globally defined fibre models. In the rest of the paper, these models are referred to as SnxRy, where S stands for Streeter-based and nxRy refers to the value of n and R used. For example, the fibre model Sn0.5R70 is obtained by imposing R~70 and n~0:5 in eqn. (1). Figure 2(a) shows the transmural variation of a for a subset of these fibre models.
The global fibre models are compared with a subject-specific fibre model developed by Karadag et al. [5]. This is based on a local optimisation of R and n carried out on 6 subjects from the DT-MRI rat dataset previously described. The cine-MRI rat scan used to generate the mesh comes from one of the subjects from this set. The left ventricle was divided into 16 regions, an adaptation of the 17-region classification of the American Heart Association [20], and on each of them eqn. (1) was fitted in order to minimise the angular difference with the data. The result is a 16-regions model of fibres angles, with a single couple (R,n) obtained for each region. We refer to this model as Karadag. Models of sheet orientation. The standard approach to the modelling of sheet orientation is that of a smooth variation in the transmural direction [2]. Recent experiments question this modelling approach showing the presence in the tissue of sheet populations with opposing orientation angles [3] and a strong level of discontinuity [13]. In order to compare the effect of the two different approaches, we have compared a simple linear transmural variation of the sheet angle b with a bimodal model. The linear model for b goes from B~{80 at endocardium to B~z80 at epicardium, and is described by b~B(2d{1). The bimodal model consists of a continuous approximation obtained by means of the sigmoidal function in eqn. ( where K is a parameter affecting the steepness of the sigmoidal function. We chose K~20 in order to have a distinct but smooth step in the trace. Figure 2(b) shows the transmural variation of the sheet angle b for the two models used in this study. In the following, unless it is specified, the linear model for the sheet angle is used.

Electromechanical model
Model of electrophysiology. At the tissue level the monodomain equations (eqn. (3) and (4)) for electrical wave propagation are solved, where V is the transmembrane voltage (mV ), t is time (ms), C m is the membrane capacitance (mF =cm 2 ), x is the surface area-tovolume ratio (cm {1 ), s is the effective conductivity (mS=cm), a matrix encoding the conductivity along the fibre, sheet and sheetnormal directions. The function I ion (g,V ) is the ionic current provided by the cell model (mA=cm 2 ) and I stim is the external stimulus (mA=cm 3 ). The vector g represents the set of state variables forming the cell electrophysiological model, whose temporal variation is regulated by the ordinary differential equation (ODE) system in eqn. (4). The model of cellular electrophysiology is that of Wang and Sobie [21]. Table 2 provides the parameters values we used for eqn. (3). The value of the surface area-to-volume ratio x is typical for rat myocytes [22]. The conductivities were adjusted in order to ensure that the resulting conduction velocity in the fibre and cross-fibre directions was within experimental range, with ratio cv fibre cv cross{fibre &2:3 [23].
Cubic wedge simulations of electrophysiology with straight fibres parallel to one of the coordinates axes were used for this purpose. Model of mechanical behaviour. The mechanical component is coupled to the electrical component via the transmembrane voltage. At the tissue level, the mechanical behaviour is modelled according to finite elasticity theory. Let us consider two systems of Cartesian coordinates, a reference system defined on the undeformed body, X~(X 1 ,X 2 ,X 3 ), and a second system representing the deformed state, x~(x 1 ,x 2 ,x 3 ). The components of the deformation gradient tensor F, which describes the transformation from the undeformed state to deformed one, can be written as C MN~P i F iM F iN , describes the deformation in the undeformed coordinate system. Finally, the Lagrangian strain tensor is The 2 nd Piola-Kirchoff stress tensor T for the total stress produced by the tissue is the sum T~T passive zT active of the stress due to the tissue passive response, T passive , and active contraction, T active . The constitutive relation between passive stress and strain is given by T passive~L W LE , where W~W (E) is the exponential strain-energy function used by Usyk et al. [4] and given in eqn. (5), where J~det(F ), C a~0 :88 kPa is a stiffness constant and b ij for fi,jg~ff ,s,ng are stiffness parameters. The parameter C b regulates the level of compressibility allowed by the material law: the higher the value the closer it is to being incompressible. After a preliminary analysis of the effect C b on the overall compressibility of the tissue we chose C b~2 00kPa, since it ensures a negligible volume reduction VR~(V end {V start )=V start ½ 100%~{1:63% between beginning and end of contraction. Furthermore the distribution of the Jacobian of the deformation J over the mesh elements between beginning and end of contraction has an average 5 th -95 th percentile range of 0:95-1:00, confirming that the model is nearly incompressible.
Two models of material passive behaviour were tested representing transversely isotropic and orthotropic behaviour. The stiffness parameters for the orthotropic case are taken from [4]. The parameters for the transversely isotropic law are obtained by averaging the relevant parameters of the orthotropic law. The parameter values for both laws are reported in Table 3.
The active stress at the tissue level, T active , depends on the active tension generated at the cell level. The cell contraction model we used was originally developed by Kerkchoffs [24] for canine myocytes. The model parameters associated with the mechanical activation time in the twitch component of the model, (b,t r and t d ), were reduced by 50% in order to get a shorter contraction duration (&150ms), more compatible with the rat cardiac cycle frequency, while still generating an analogous amount of active tension (&16kPa at peak contraction) in that shorter interval. Finally, the total stress T is required to satisfy the conservation of momentum as shown in eq. (6) where the effect of inertia and gravity is neglected. The system is solved with appropriate Dirichlet boundary conditions for displacement and Neumann boundary conditions for tractions, in this case the blood cavity pressure onto the endocardial surface, which is defined as follows.
Cavity pressure trace. The temporal evolution of the intracavital blood pressure is approximated by means of a sigmoidal function ranging from {1 kPa to {5 kPa, where the negative sign is due to the fact the pressure is applied in the inward direction along the normal to the surface. The range of pressure applied is significantly lower than that typically measured in the rat left ventricle, which reaches approximately 16 kPa (or more, for example during physical activity) at end-systole [25]. This is a known limitation of some contraction models and follows existing computational modelling literature [26]. However, the chosen pressure trace ensures that the left ventricle undergoes first an isovolumic contraction and then an ejection phase. Figure 3(a) shows the pressure trace with an example of the resulting cavity volume.

Electromechanical simulations
All simulations were run using Chaste, an open-source simulation package for computational biology developed in our group (http://www.cs.ox.ac.uk/chaste) [27]. The simulation duration is 150ms which is approximately the time at which full contraction is achieved. The timesteps used for solving the electrophysiology are 5 : 10 {3 ms for the cell model (eq. (4)) and 1 : 10 {2 ms for the monodomain equation, while for the mechanics they are 0:1 ms for the contraction model and 1 ms for eq. (6). The numerical scheme used has been described in previous work [28].  Table 3. Stiffness parameters values for the material law (eq. (5)).

Global and local metrics
The global measures used to evaluate the simulations are: ejection fraction EF , wall thickening WT and shortening in the base-to-apex direction SL. All the quantities are computed by comparing the state of the mesh at the beginning (t 0~0 ms) and end (t f~1 50ms) of contraction. EF is calculated as the relative difference between the mesh cavity volume vol at t 0 and t f , , where w at a given time is obtained by averaging four measures of wall-thickness taken symmetrically onto a short-axis slice in the midventricular region (at approximately 25% of the long axis from the base). We chose this specific short-axis slice to compare with the cine images, where WT is more easily measured in the midventricular region. SL is calculated by comparing the z coordinate of the epicardial apical node along the z axis (which coincides with the long axis of the left ventricle, see Figure 1(a)) at t 0 and t f , . This is an adequate measure of the shortening along the base-to-apex direction since the base is kept fixed in the model and the apex moves almost exclusively along the long axis (of course, in situ, the movement is in the opposite direction, i.e. the atrio-ventricular valve plane moves 'downwards' in the direction of the apex, but in terms of changes in apico-basal length, this parameter offers a satisfactory approximation within the modelling framework used). Here SLv0 represents shortening, while SLw0 represents elongation of the ventricle during contraction.
In order to compare the predicted values with experimental data, we measured these metrics on the cine-MRI dataset used for building the ellipsoidal mesh by applying analogous methods. The measured values are EF cine~7 7%, WT cine~5 4% which are consistent with values in the literature [29], and SL cine~{ 10%.
As local metrics we considered the Lagrangian circumferential strain Ecc, and two shear strains, the circumferential-longitudinal shear strain Ecl, and the circumferential-radial shear strain Ecr. The Lagrangian strain tensor E in the Cartesian system is obtained from the deformation tensor F calculated during the simulation, that was sampled every 10 ms. Then the conversion to cylindrical coordinates was applied. We could not compare these strain predictions directly with data, but they show good agreement with results in the literature as further explained in the Discussion.

Results
In this section we present the results of the series of electromechanical simulations we performed in order to (a) compare sixteen globally defined fibre models (SnxRy for x~0:5,1,2,3 and y~50,60,70,80) against one subject-specific fibre model (Karadag), and (b) evaluate the effect of introducing sheets by comparing transverse isotropy to orthotropy and with linear and bimodal sheet models. The results are grouped into global and local measures.

Global measures
The global measures considered in this study are 1) ejection fraction EF , 2) wall thickening WT, and 3) shortening in the longaxis direction SL. The variation in EF due to the different fibre models is shown in Figure 3(b), with the transversely isotropic material law. EF is plotted vs R, and the global fibre models SnxRy are grouped according to the exponent n in order to highlight the pattern of EF with increasing values of R for each value of the exponent n~0:5,1,2,3. The single value of EF for the Karadag fibres is represented as an horizontal line. The predicted EF is overall in the range 50%-60%, which can be considered physiological [30], but which is still lower than the measured value EF cine~7 7%. Among the global fibre models, for n~1,2,3 increasing values of the angle R and decreasing values of the exponent n raise the value of EF . The reverse holds for n~0:5, for which the increase in the maximum angle R causes a rapid decrease of the EF . A possible explanation for this downward Changing from the transversely isotropic to the orthotropic material law (with the linear sheet model) has the general effect of increasing the EF value for all the fibre models, with an average increase of 2:7 pp (percentage points), resulting in EF values in the range 53%-63%, i.e. still below the measured value. Figure 4(a) shows the comparison of EF between transversely isotropic and orthotropic material laws with a linear sheet model for the simulations with global fibres with n~1 and Karadag fibres. Similar results are obtained for the other values of n. Within the orthotropic case, the effect of using a linear or a bimodal sheet model has the interesting effect of producing predictions much closer to those obtained under transverse isotropy than to orthotropy with the linear sheet model, but in general with a less consistent pattern across the fibre models (and therefore this is not shown).
For WT, Figure 5(a) represents the wall thickening WT vs R, with the transversely isotropic material law. For the global fibre models, we observe the same trend as with EF , i.e. for increasing values of R and decreasing values of n WT increases. Karadag fibres still produce one of the highest values of WT, but some global fibres models also give good results. The use of the orthotropic material law with the linear sheet model further increases the value of WT for every fibre model, as shown in Figure 4(b), with an average increase of 4:3 pp. The effect of varying the sheet model from linear to bimodal produces again values closer to the transversely isotropic case. Overall the range of variability for WT is 16%-37%, with only some fibre models giving values close to the physiological range. In any case there is a clear discrepancy between the predicted WT values and that measured on cine-MRI, WT cine~5 4%.
The shortening in the base-to-apex direction (negative values of SL) is also substantially affected by the fibre type and material law, ranging from {15% to z10%. Figure 5(b) shows SL vs R, with the transversely isotropic material law. Within the global fibre models SnxRy, high values of the exponent n combined with low values of R produce lengthening (positive values of SL), whereas there is increasing shortening for low values of n and increasing values of R. In particular, some global fibre models match or exceed the measured shortening SL cine~{ 10%. Karadag fibres produce a shortening value smaller than the measured one. The effect of changing form transversely isotropic to orthotropic material law (with the linear sheet model) is shown in Figure 4(c).
The orthotropic material law increases the shortening by 3:5 pp on average. By changing the sheet model to bimodal approximately 1 pp less shortening was obtained.

Local measures
We considered the circumferential strain Ecc, the circumferential-longitudinal shear strain Ecl and the circumferential-radial shear strain Ecr, which were uniformly averaged in the basal (upper 35% of the left ventricle), midmyocardial (middle 35% of the left ventricle) and apical sections. We show only the average strain in the midmyocardial region, because it is a good measure of the average behaviour of the left ventricle. In the absence of direct strain data, we compared the predicted traces with experimental results in the literature. Figure 6(a) shows the Ecc trace for the subset of global fibres models with n~1 and for the Karadag model in the case of transverse isotropy. The Karadag model has an intermediate behaviour between the fibre models Sn1R70 and Sn1R80, as we have already found for the global measures. Amongst the global fibre models, the magnitude of the strain decreases for increasing values of the angle R. All traces start from a positive value of Ecc rather than zero, due to the passive expansion performed prior to the electromechanical simulation, which leaves a residual strain component in Ecc. In order to compare predicted Ecc traces with experimental results on Ecc in rat, the study by Daire et al. [29], who employed cine-MRI and Tagged-MRI, was considered. The measurements of Ecc were taken by Daire et al. [29] at fixed percentages of the cardiac cycle, therefore it was necessary to convert this temporal subdivision into a sequence of time points in the range 0-150ms (the duration of our simulations). The data (mean and standard deviation) for the first half of the cardiac cycle, which corresponds approximately to the simulation duration, are shown in red in Figure 6(b). The green line is the Ecc trace for the simulation with the Karadag fibre model and transverse isotropy, after subtraction of the residual strain from the pre-simulation passive expansion. Apart from the inevitable inaccuracy of the conversion from cardiac cycle percentages to the simulation time interval, there is very good agreement between predicted and measured Ecc.  heart was not available in the literature. Instead, we compared predicted shear strains with the measures performed in mouse by Zhong et al. [31]. As with the Ecc, it was necessary to convert from percentages of cardiac cycle to time points in the range 0-150ms. Figures 7(b) and 8(b) show the comparison between data and predicted Ecr and Ecl, respectively, for the simulation with the Karadag fibre model and transverse isotropy. The agreement for Ecr is poorer than for Ecl. However, in both cases the trends are similar to the experimental ones.
The effect of sheets was negligible on Ecc, while it was substantial on the shear strains. Changing from transverse isotropy to orthotropy, and within the orthotropic behaviour, from a linear to a bimodal model for the sheets, caused a consistent increase in the magnitude of Ecr and Ecl, as shown in Figure 9(a) -9(b) for the Karadag fibre model in the midmyocardial region.

Discussion
In the previous section, results are presented of a systematic computational study into the effect of tissue microstructure on contraction, using a simplified model of rat left ventricle. It consists of a series of electromechanical simulations carried out on a thickwalled truncated ellipsoid, with the aim of quantifying the combined effect of different models for fibre and sheet directions with different tissue passive behaviours during simulated active contraction. Sixteen global fibre models were compared against a subject-specific model. The impact of sheets was also investigated by comparing transversely isotropic and orthotropic material properties as well as two models for the sheet orientation, namely a linear and a bimodal model.
In the development of the model we tried to be as speciesspecific as possible in all the components, with the exception of the  reduced range for the cavity pressure, as described in Methods and further discussed later on. The geometry has dimensions taken from rat cine-MRI scans, the parameters of tissue electrophysiology are based on measurements in rat. Regarding the AP model, we preferred the mouse model of Wang and Sobie [21] over other available rat models for reasons of computational efficiency. Of course, the ventricular activation sequence with simulated apical pacing differs from normal in vivo patterns that were present during cine-MRI. However, as we compared pre-activation (diastole) with peak contraction (systole), we believe that this is unlikely to have had a major effect on our findings related to steady-state comparisons. As for the mechanics, there are rat-specific biophysical contraction models. However, we preferred to use the model of Kerckhoffs et al. [24] as it allowed us to adjust its parameters in a controlled manner while retaining a good description of the contractile behaviour of the myocyte. The stiffness parameters of the exponential constitutive law (eqn. (5)) for material orthotropy are taken from a study carried out on canine hearts by Usyk et al. [4]. In their paper Usyk et al. explore the stiffness parameters space in order to study their effect on mechanics, providing a range of variation for each of the parameters. We chose the orthotropic parameter set that showed magnitudes comparable with the values reported for rat in transverse isotropy by Omens et al. [32]. This was necessary because we did not find published orthotropic stiffness parameters values for exponential constitutive laws for rat. We obtained the stiffness parameters for transverse isotropy by averaging the orthotropic ones in the sheet and normal directions. We ensured  Effect of Tissue Microstructure on Contraction PLOS ONE | www.plosone.org that the two parameters sets for orthotropy and transverse isotropy had the same magnitudes in order to minimise the effect of parameters magnitude and explore only that of material symmetry.
Various measures were considered, grouped into global and local parameters, in order to study the effect of the fibre model, material law and sheet model on contraction. For the global measures, ejection fraction (EF ), wall thickening (WT) and baseto-apex shortening (SL) were compared. We chose these measures because they are representative of the main features of cardiac function in the left ventricle. Furthermore, most electromechanical models do not currently predict WT and SL well enough. The global measures are consistently affected by the fibre model and material law, while they show smaller variation when changing the sheet model. The only exception seems to be the downward trend of EF in the case of n~0:5 (Figure 3(b)). We believe this is due to an inverse correlation between the magnitude of SL and the rate of increase of EF . In other words, the more SL increases in magnitude for a given value of n and increasing values of R, the smaller the increase in EF . This effect can be seen in Figure 3(b), where the steepness of the EF curves decreases consistently going from n = 3 to n = 0.5. This observation suggests that there is a trade-off between shortening and efficiency of the pumping function. Too much shortening is as detrimental to pumping function as elongation. Regarding the fibre models, the subjectspecific fibre model does not produce values of the global measures that are notably closer to the measured values than those obtained using the best global fibre models. In particular, among the fibre models, those with low exponent n and high torsion angle R (such as Sn1R70 and Sn1R80) produce values of the global metrics that are comparable or better than those of the subject specific model. These results are encouraging in that they confirm that, in the case of the Streeter-like fibre models (see eqn. (1)), the exponent n~1 describes fibres in rat left ventricle better than other values, in accordance with previous comparisons of rule-based versus DT-MRI-derived fibres in rat [33]. In general, the predicted values for EF and WT, although being at the physiological range ( [30], [13]), are lower than the measured ones, while some global fibre models are able to match or exceed the measured value for SL. Changing from transversely isotropic to orthotropic material properties (with the linear sheet model) increases the predicted values of the three metrics, but this is not enough to match the measured values for EF and WT. Within the orthotropic setting, no significant effect was seen on SL, while for EF and WT the bimodal model case gave results intermediate between those for transverse isotropy and orthotropy with the linear sheet model.
The Lagrangian strains Ecc, Ecl and Ecr are used as local measures. The Karadag subject-specific model has strain traces close to those of global fibre models with lower exponent n and higher angle R. In particular the models Sn1R70 and Sn1R80 are those that in terms of strains give the most similar results to the Karadag model. This confirms our findings in the analysis of the global metrics. Interestingly, the fibre models that produce higher magnitude of strains (that is higher n and lower R) are the same that gave poorest results in terms of the global metrics. Conversely, this means that fibre models such as Karadag, Sn1R70 and Sn1R80 are able to convert smaller strains into more efficient cardiac function. Changing from transversely isotropic to orthotropic material properties has a negligible effect on Ecc, a small effect on Ecl, while a more pronounced change was evident in Ecr. In particular, introducing the bimodal model increases the average magnitude of Ecr in the midmyocardial region by w30% with respect to both transverse isotropy and orthotropy with the linear sheet model. This is as expected since the effect of introducing orthotropy is more likely to have a larger effect on Ecr rather than Ecc. In general, the results show strain traces and magnitudes close to those reported in the literature for rodents [29] [31], especially in the case of Ecc, as shown in Figures 6(b), 7(b) and 8(b), for Ecc, Ecr and Ecl traces, respectively, for the simulation with Karadag fibre model and transverse isotropy.
It is important to discuss the limitations of our model. (a) The use of a simplified geometry: the natural geometrical variability of the left ventricle is not taken into account in this study, including the lack of a neighbouring contracting chamber (right ventricle), or of papillary muscles whose orientation is predominantly axial. Previous studies suggest this might not bring a significant improvement [9], nonetheless, it would be interesting to investigate if the conclusions about the subject-specific vs global fibre models hold also in the case of a subject-specific geometry. (b) The electrical activation pattern: the apical region was chosen for the electrical stimulation in order to mimic the apex-to-base direction of the electrical activation in the left ventricle. This type of activation pattern is consistent, for example, with the type of external stimulation employed in the Langendorff-perfused experiments. However, this is not what happens in the living heart [34]. Evaluating the influence of changes in activation pattern on predicted left ventricular mechanics would be an interesting application of our model. (c) The ad-hoc procedure to estimate the pre-P mesh prior to applying the pressure load: this step should be tackled in a more automatic and rigorous approach, by solving an inverse problem. (d) The reduced range of pressures applied in the model: this is probably the most important limitation of the model as pressure is one of the best characterised physiological parameters of the heart. If we compare the results of the preliminary study we carried out using the model without pressure [17] with those presented here, we find similar trends in the global measures for the different fibre models. By introducing pressure in the model, we find that the magnitude of those measures increases significantly on average, meaning that the pressure boundary condition is fundamental for the efficiency of the simulated pumping function, but does not interfere with the effect of the different fibre models on contraction. The use of a reduced pressure range is the result of the best compromise between realistic contraction patterns and physiological pressures in the model. With the current modelling setting, the use of realistic (based on laboratory measurements) material properties in combination with the chosen contraction model forces the pressure to be reduced to avoid non-physiological contraction patterns. The alternative would be to enforce physiological pressures by modifying the material properties (which would no longer be based on experimental measurements) and fully reparametrizing the contraction model (for which experimental data would not be available). We chose the first approach for several reasons. The first is to reduce the degrees of freedom of the model, setting one parameter (the pressure) rather than a set of mechanical properties. Another is that of consistency with existing literature, which allows us to compare results directly. Still, as mentioned above, this is a major study limitation. As such, it identifies a key direction for further study, reinforcing an important role of quantitative models as a drivers of scientific progress [35], and highlighting that currently available data, implemented mechanisms, and attempted explanations are insufficient to fully reproduce reality. (e) The lack of the transverse angle in the fibre model definitions: this angle regulates the level of crossover of fibres which may be important in a correct prediction of shear strains [16]. (f) Cross-fibre generation of active tension: at the tissue level, the model generates active tension along the fibre direction only. Experimental studies have shown that there is a certain amount of active tension generated in the cross-fibre plane as well [36], and some important electromechanical simulation studies already take this into consideration [4], [16]. In particular, Bovendeerd et al. [16] have shown that including the cross-fibre active tension generation into their model improved the prediction of Ecr. (g) The fibre models considered: the choice of rule-based and subject-specific fibre models analysed was necessarily reduced in order to focus on a specific scientific question. It is important to note that we did not consider asymmetrical fibre angle intervals for the rule-based model, which have been found experimentally [1] and might affect, at least locally, the predicted contraction patterns. However, we expect that their effect would remain within the ranges predicted by the fibre models already examined.
Eriksson et al [37], for example, have compared symmetric and asymmetric transmural fibre angles for a linear rule-based fibre model and found generally consistent predictions in the metrics considered. Another important aspect we did not take into account is the variability in subject-specific fibre distribution. Karadag et al [5] report standard deviations of + 15 degrees for R and + 0.25 for n. It would be interesting to investigate the inter-subject variability of the Karadag fibre model and its effect on contraction patterns in the future. (h) The sheet orientation models: we used two simplified models of sheet transmural variation, neglecting the great regional variability in terms of sheet patterns that has been found experimentally, and that is thought to give rise to a concertina-like re-arrangement during contraction and relaxation [13]. Overall the impact of sheet inclusion in our study was small, with the exception of changes in Ecr. It could be argued that this could be due to the high stiffness in the sheet direction imposed by the constitutive law. However, preliminary tests with different magnitudes of the stiffness parameters did not show such a substantial variation in the global metrics results to make us believe that this problem could be addressed simply by tuning the constitutive law parameters. Furthermore, introducing the more realistic bimodal sheet model did not affect the global metrics in the way we expected. For example, it did not cause a significant increase in the predicted wall-thickening, while it has been shown that this mechanism is directly related to sheets dynamics [14].
To conclude, we developed and tested an electromechanical model and performed a computational study of the impact of tissue microstructure on contraction of the left ventricle in rat. We showed that the model predicts values of global metrics that are within physiological ranges but, except for SL, they fall below the values measured on cine-MRI scans. Introducing subject-specific fibre models, as opposed to global fibre models, does not provide a substantial increase in EF and WT values. The analysis of strains confirmed that the subject-specific fibre model gave predictions similar to linear global fibre models with higher rotation angle. Changing from transverse isotropy to orthotropy in the passive behaviour, therefore introducing the effect of tissue anisotropy due to the sheets, improves the predictions of EF and WT, but they still remain lower than the measured values. This means that including the sheet orientation in the model is an important factor, and a more detailed representation of angular re-orientation of sheet material may be required. Therefore, we believe that, within the current modelling framework, the subject-specific approach in fibre definition by itself is not able to compensate for the model limitations that are believed to come from the underestimation of sheet dynamics, and that further investigation of the dynamics of sheet orientation should be carried out and novel approaches developed in order to fully represent its effect on cardiac mechanics.