## Abstract

A scientific understanding of individual variation is key to personalized medicine, integrating genotypic and phenotypic information via computational physiology. Genetic effects are often context-dependent, differing between genetic backgrounds or physiological states such as disease. Here, we analyse *in silico* genotype–phenotype maps (GP map) for a soft-tissue mechanics model of the passive inflation phase of the heartbeat, contrasting the effects of microstructural and other low-level parameters assumed to be genetically influenced, under normal, concentrically hypertrophic and eccentrically hypertrophic geometries. For a large number of parameter scenarios, representing mock genetic variation in low-level parameters, we computed phenotypes describing the deformation of the heart during inflation. The GP map was characterized by variance decompositions for each phenotype with respect to each parameter. As hypothesized, the concentric geometry allowed more low-level parameters to contribute to variation in shape phenotypes. In addition, the relative importance of overall stiffness and fibre stiffness differed between geometries. Otherwise, the GP map was largely similar for the different heart geometries, with little genetic interaction between the parameters included in this study. We argue that personalized medicine can benefit from a combination of causally cohesive genotype–phenotype modelling, and strategic phenotyping that captures effect modifiers not explicitly included in the mechanistic model.

## 1. Introduction

The emerging science of personalized medicine will have to be founded on computational physiology and a scientific understanding of individual variation [1]. This variation results from a complex interplay of genetic, age-related, environmental and lifestyle factors. Its unravelling requires the incorporation of computational physiology into genetics theory, a unification that has only just begun [2–5].

In genetics, a *phenotype* is any observable trait of interest in an organism, whereas its *genotype* denotes its hereditary material or a relevant subset thereof. The relation between genotype and phenotype can be conceptualized as a *genotype–phenotype map* (GP map), assigning a phenotype to each possible genotype conditional on the environment [5]. However, standard population genetic models simply assign phenotypic values directly to genotypes, without any intervening causal explanation. In contrast, computational physiology is capable of mimicking variation in measurable phenotypes based on mechanistic understanding of dynamical systems [6]. It is a small conceptual step to map genotypes to low-level parameters of physiological models, rather than directly to the high-level phenotypes that are the eventual target of understanding [4,5].

The basic premise of *causally cohesive genotype–phenotype* (*cGP*) modelling [2,5] is that in a model that is capable of accounting for the phenotypic variation in a population, the causative genetic variation will manifest in the model parameters. Thus, a cGP model denotes a model where parameters at some level have an articulated relationship to the individual's genotype, and higher-level phenotypes emerge from the mathematical model describing the causal dynamic between the lower-level processes. In this view, individual variation corresponds to different points in parameter space (for quantities that are constant over the timescale being studied) and in state space (for quantities whose dynamics are being explicitly modelled). Furthermore, what is treated as a parameter in one model can usually be conceived as the output of more fundamental mechanistic models [2].

The mammalian heart is an excellent cGP model system [3], being modelled at levels from ion channels to the whole organ. In a previous study of a single-cell model [3], we demonstrated the deep analogy between many classical genetic concepts and sensitivity analysis of cGP models, describing the sensitivity of high-level and intermediary phenotypes with respect to genotype, via genetically influenced model parameters. Linear main effects in a regression model are akin to intra- and interlocus additive gene effects; nonlinear effects are akin to genetic intralocus dominance; interaction terms are akin to interlocus epistatic interactions (i.e. where the effect of a genotypic substitution at one locus depends on the genotype at another locus). Such epistatic interactions or genotype × environment interactions can result in incomplete penetrance for binary phenotypes such as disease, where not all carriers of a predisposing allele may actually express the disease. The research field of sensitivity analysis is well established, with a large mathematical and practical literature [7–11], and is thus an important operational tool for linking genetic concepts to system dynamics models. For example, we recently showed that *in silico* genome-wide association studies (GWAS) give results similar to a sensitivity analysis of a cGP model [4]. Thus, a sensitivity analysis can guide GWAS analyses by indicating where genetic variation is more likely to propagate to specific phenotypes.

Multiscale heart models offer additional challenges than single-cell models, owing to the added levels of organization (tissue, organ), spatiality of processes and phenotypes and branches of physics. Many genetically influenced cell-level properties have been linked to biomedically important phenomena such as morphogenesis, remodelling and hypertrophic growth (table 1 and [15]). As reviewed in [16], continuum mechanics is a useful framework for bridging the genotype–phenotype gap in cardiac disease. Furthermore, features of the GP map may change in disease [17,18], as the organism is brought outside its most robust domain of normal operation.

Here, we investigate the extent to which GP map features of the heart as a soft-tissue mechanics system are modified by pathological changes in heart geometry, namely concentric and eccentric hypertrophy (figure 1). Specifically, we model the passive filling phase (late diastole) of the left ventricle, in which many individual factors have been studied previously [22–25]. The mechanics in this phase are relatively simple owing to the absence of active contraction, yet the strain and elastic energy stored in diastole gives our conclusions some relevance to the later phases of the heartbeat as well. Owing to the simple mechanics in the filling phase, the GP map is likely to be simpler than for a full heart cycle with coupled electromechanics. However, we hypothesize that more of the low-level parameters (e.g. fibre angles and stiffness of cardiac tissue) will influence the higher-level phenotypic outcome (e.g. change in shape and volume of the heart) for the concentric geometry, because concentric hypertrophy primarily impairs filling, whereas eccentric hypertrophy impairs contraction. Thus, this study emphasizes a phase of the heartbeat whose problems are more relevant to concentric than eccentric hypertrophy.

## 2. Models and methods

### 2.1. Mechanics

The framework for mechanics is described in Land *et al.* [26]. We ran simulations of the passive inflation of the ventricle where the ventricular deformation satisfies the force balance of the myocardial elastic response as a function of left ventricular pressure (LVP). Before deformation, i.e. at zero pressure, the inner and outer surfaces were set to be truncated ellipsoids (figure 1), as done in several other studies [23,27–29].

The myocardium was modelled as a hyperelastic, orthotropic and incompressible material [30]. We formulate the stress equilibrium equations by using the principle of virtual work on a general deformable body with initial volume *V* and surface area *A*
2.1

Here, **S** is the second Piola–Kirchhoff stress tensor and *δ***E** is the virtual Lagrangian strain tensor. *p*_{app} is the applied ventricular pressure acting on the inner wall (endocardial) surface, **F**^{–T} is the inverse of the transposed deformation gradient tensor, **N** is the normal vector for the undeformed material surface, whereas *δ***u** is an arbitrary virtual displacement. For hyperelastic materials, the components of the second Piola–Kirchhoff stress tensor are calculated by differentiating the strain energy function by the corresponding components of the Lagrangian strain tensor. For incompressible materials, the strains are separated into distortional and dilatational components, and we calculate the second Piola–Kirchhoff stress tensor from
2.2where is the strain energy function and **F** is the deformation gradient tensor. The scalar *p* serves as a Lagrangian multiplier interpreted as the hydrostatic pressure in the material/tissue. The work done by the hydrostatic forces is zero owing to the incompressibility constraint, expressed on variational form as
2.3where *J* is the determinant of the deformation gradient tensor **F**. The stress–strain relation in hyperelastic models follows from our assumption of a strain energy function, . The Costa law [31] describes an orthotropic tissue with different stiffness in the three orthogonal directions (see §2.3 Microstructure),
2.4Here, *a* is an overall stress scaling parameter, and *Q* is given by
2.5

In this relation, *E _{ij}* and

*b*are respectively the Lagrangian strain components and stiffness parameters associated with the microstructural directions

_{ij}*s*,

*f*,

*n*(see Microstructure section).

### 2.2. Geometry

The geometries are shown in figure 1, and corresponding parameter values are presented in table 2. The hypertrophic and dilated geometries capture the essential features of observed pathological deviations from normal heart shapes [32,33] within the simple model of a truncated ellipsoid. We note that assuming such a simplified heart geometry means that the simulated deformation is probably more representative of the free wall than of the septum.

### 2.3. Microstructure

The microstructure of the myocardium is characterized by the fibre (*f*), sheet (*s*) and sheet-normal (*n*) directions. The microstructure of any material point in the myocardium is identified by means of three numbers; the fibre angle, the sheet-angle and the imbrication angle. We use the same definition of fibre, sheet and imbrication angles as in [34]. Here, the fibre angle was set to vary linearly transmurally between the prescribed input angles *α*_{endo} and *α*_{epi} at the inner and outer wall, respectively. We incorporate sheet structure based on the experimental findings reported in [35]. In order to describe the sheet angle we used a single parameter *β*, where . The sheet angle was linearly interpolated between these values. The imbrication angle was set to zero throughout the tissue for all simulations. Figure 2 shows an example on how the fibres and sheets angles vary throughout the ventricular wall in the model.

### 2.4. Simulation of passive filling

The simulations of the passive inflation of the left ventricle from 0 to 10 mmHg (1.33 kPa) in 11 equally spaced pressure steps were performed in a cubic Hermite interpolated finite-element model with 8 × 7 × 2 elements. Only parameter combinations that converged for all applied pressures were analysed.

### 2.5. Boundary conditions

We defined the following boundary conditions: within the left ventricle, we imposed a uniform LVP by inserting a non-zero normal traction vector into equation (2.1) along the inner wall of the inner elements. At the base of the ventricle, the nodes were fixed, so no basal material points were allowed to move out of the plane, *z* = 0. In addition, the endocardial basal nodes which were aligned with the *x-* and *y*-axis in the reference configuration were only allowed to move radially.

### 2.6. *In silico* genetic variation

Although many important parameters and phenotypes for heart mechanics have been shown to have a genetic basis (table 1), the details of that basis remain unknown. We therefore make very simplifying assumptions about the gene-to-parameter map, as in [3], and also assume that geometry and genetic parameters are independent. For each geometry, we introduced genetic variation by varying 10 parameters of the model; the seven stiffness parameters of the Costa law (*a*, *b _{ff}*,

*b*,

_{ss}*b*,

_{nn}*b*,

_{fs}*b*,

_{fn}*b*) and the three microstructural parameters for the fibre orientation (

_{sn}*α*

_{endo},

*α*

_{epi}and

*β*). As in [3], we assumed that each parameter was determined by one locus with ‘low’ and ‘high’ alleles A and a, with genotypic values of aa = 50%, Aa = 100% and AA = 150% of the baseline parameter estimate (table 3). Baseline values for the stiffness parameters were obtained from [36].

We performed a full factorial design for every combination of four of the 10 parameters in table 3 (keeping all other parameters at their base value; giving 4521 parameter sets). The rationale for this design was that it allowed us to explore up to fourth-order interactions between parameters that describe the stiffness and microstructure.

### 2.7. Phenotypic measures

We focused on two main groups of organ-level phenotypic traits associated with the passive inflation of the left ventricle, the first group characterizing the end-diastolic deformed geometry of the ventricle (shape phenotypes) and the second group characterizing the stress– and strain–energy fields within the myocardium (field phenotypes). The shape phenotypes included the change in ventricular volume (Δ*V*), the change in long- and short-axis diameters (Δ*La* and Δ*Sa*), the mechanical energy stored in the current configuration given by the area below the pressure–volume curve , the change in wall thickness (Δ*WT*), calculated from the volume of the tissue divided by the average area of the inner and outer wall, and the inner and outer twist *φ*_{in} and *φ*_{out}. The inner and outer twist was defined as the difference of angular displacement (measured in radians) between the uncollapsed (endo- and epicardial) nodes of the apical elements and the basal (endo- and epicardial) nodes.

Furthermore, incorporating spatial mechanics requires phenotypic measures that summarize the spatial distribution of the stress and strain fields embedded in biological tissue. These field phenotypes included the spatial mean and variance (calculated as volume integrals) of the Cauchy fibre stress (*σ _{ff}*), the Langrangian fibre strain (

*E*), the strain energy () as well as the proportion of

_{ff}*Q*attributable to the fibre direction: . The strain energy in the Costa law is a function of all the different strains and the stiffness parameters (unlike e.g. in the pole-zero law [37]). Hence, it is not possible to decompose the total stored elastic energy into the different microstructural directions, and we therefore use

*Q*as a measure of the proportion of elastic energy stored along the fibres.

_{ff}### 2.8. Characterizing the genotype–phenotype map

Based on the end-diastolic deformations, we used multiple linear regression where the stiffness and microstructure parameters were independent continuous variables, forming an (*n* = 4521 by *K* = 10) matrix **P**, and the phenotypes were dependent variables. In other words, for each phenotype *φ*, we fitted the model
2.6where (*p _{i}*

_{1},

*p*

_{i}_{2}, …) is the parameter scenario with index

*i*, and

*ɛ*is the residual, i.e. the difference between the phenotypic value from the physics simulation and the statistical approximation.

_{i}Because the columns of **P** are orthogonal to each other (by virtue of the experimental design) as well as to the residual vector ** ɛ** (by virtue of least-squares fitting), the variance in

*φ*can be decomposed into residual variance plus variance components for each genetic parameter.

In the widest sense, we use an *interaction* between the *k*th and *l*th parameter (in their effect on *φ*) to mean that the effect of a change in the *k*th parameter depends on the value of the *l*th parameter and vice versa. One could thus also call the *l*th parameter an *effect modifier* of the *k*th parameter. Such interactions can be represented by forming new terms by multiplying the original predictors, e.g. adding *β*_{kl}*p _{ik}p_{il}* to the regression equation. This effectively replaces

*β*with (

_{k}*β*+

_{k}*β*)

_{kl}*p*or, equivalently,

_{il}*β*with (

_{l}*β*+

_{l}*β*)

_{kl}*p*, thus making the effect of one

_{ik}*p*depend on the value of another.

Nonlinear effects can be represented by treating the *k*th parameter as a factor with levels (*aa*, *Aa*, *AA*), e.g. using *Aa* as the reference level and replacing *β*_{k}*p _{ik}* by

*β*

_{aa}

*D*+

_{aa}*β*

_{AA}

*D*, where the dummy variables (

_{AA}*D*,

_{aa}*D*) are (1, 0), (0, 0) and (0, 1) for genotypes

_{AA}*aa*,

*Aa*and

*AA*, respectively.

The proportion of phenotypic variance explained by each genetic parameter indicates its main effect in the GP map. Large residual variation would indicate phenotypes whose GP map included strong nonlinearities or parameter interactions. The set of variance decompositions for all the phenotypes forms a basis for comparison of GP map features between the normal and pathological heart geometries.

Scripts and data to reproduce the figures are provided in electronic supplementary material.

## 3. Results and discussion

### 3.1. Genotype–phenotype map features for normal versus pathological heart geometries

The phenotypic effect of the simulated parameter variation is summarized in figure 3, and the variance decompositions in figure 4 show which parameters have the strongest impact on each phenotype. Together, this gives an overview of the simulated GP map for the normal and pathological heart geometries. Overall, pathological geometry did little to modify the GP map, as reflected in the variance decompositions (figure 4). Although the mean value of most phenotypes differed strongly between geometries, there was little impact on relative variability (height of boxplots in figure 3, measured on a logarithmic phenotype scale).

The GP map was strikingly linear for most shape phenotypes, as a multiple linear regression analysis accounted for more than 90% of the phenotypic variation (figure 4*b*). However, phenotypes describing twist and spatial variability in strain had larger residual variance, indicating either nonlinear ‘gene’ effects or epistatic interactions. Allowing for nonlinear effects of microstructural parameters by treating them as categorical explanatory variables (figure 4*c*) accounted for most of this, although the spatial strain variability phenotypes retained about 20% residual variance, which can be attributed to interaction.

The GP map in a cGP study is a composite of the genotype–parameter map and the parameter–phenotype map. These two submappings should be considered separately when discussing the low proportion of nonlinearity and interactions in figure 4. Uncertainty and variability [25] and interdependency [38] of material parameters are fundamental challenges in biomechanical modelling. Here, the genotype–parameter map was modelled additively, meaning that all residual variance in figure 4 results from the physiological model that maps parameters to phenotypes. Given that biological systems are inherently nonlinear, the assumption of an additive mapping from one gene to one parameter is clearly unrealistic. However, additive genotype-to-parameter maps is a natural starting point for most exploratory cGP studies [3,39,40], because then any nonlinearity that emerges must be due to the modelled physiology. But, when the biomechanical model is clearly nonlinear, why do we not observe a more nonlinear GP map? One reason may be that the nonlinearities found, for instance, in the constitutive laws in this study, are *monotone* relationships. Gjuvsland *et al.* [41] argued that systems with monotone dose–response relationships tend to produce more additive GP maps than unimodal or other non-monotone relationships. The low amount of interaction effects in figure 4 suggests that findings from earlier separate studies of stiffness [22], geometry [23] and fibre structure [24] will often remain valid under different conditions, including different genetic backgrounds and more complex and realistic genetic parameter variation.

The details of the parameter-to-phenotype mapping will depend on the choice of constitutive law, as they use different parametrizations to capture different aspects of material behaviour. However, inasmuch as many constitutive laws span roughly the same repertoire of deformation behaviours [42], we would expect general conclusions to hold about the relations of *groups* of parameters (e.g. stiffness parameters versus microstructural parameters) to certain phenotypes. For the transversally isotropic Guccione law [43], we would expect quite similar results as in this study, because the Costa law has the Guccione law as a special case, and most phenotypes were rather insensitive to the parameters that distinguish the Costa law from the Guccione law.

### 3.2. Propagation of variability

In agreement with our hypothesis, the concentric geometry did allow more parameters to propagate their variation to organ phenotypes (‘shape phenotypes’ in figure 4*c*, middle). The main difference between the concentric and eccentric scenarios was the relative importance of the fibre stiffness parameter, *b _{ff}* and the overall stiffness parameter,

*a*(figure 4). Fibre stiffness was a chief explanatory variable for most phenotypes except twist and fibre stress (

*σ*), but relatively less important in the concentric case. Overall stiffness was important for the shape phenotypes, but less so for the eccentric case. A likely explanation is the occurrence of the term

_{ff}*b*

_{ff}E^{2}

*in the exponent of the constitutive law (equation (2.5)). Strains are generally smaller for the concentric geometry (figure 3), and so there is less effect of variation in*

_{ff}*b*in this case. This can be given a protein-level interpretation owing to the contribution to passive tension from titin versus collagen. Tension is chiefly owing to titin at small strains up to about sarcomere length of 2.0 μm, whereas at a sarcomere length of 2.2 μm, titin and collagen contribute equally [44]. This suggests that the genetic and physiological basis of

_{ff}*a*may involve titin, whereas that of

*b*may be more owing to collagen. Furthermore, concentrically hypertrophic hearts may be more sensitive to variation in the expression and properties of titin, whereas collagen may be more important to eccentrically hypertrophic hearts.

_{ff}There was a striking lack of propagation to whole-organ phenotypes from the variation in sheet structure (*β*) (in agreement with Wang *et al*. [24]) and in most individual parameters of the Costa material law (*b _{fs}*,

*b*, … ; figure 4). The insignificance of

_{fn}*β*runs contrary to the conclusions of LeGrice

*et al.*[45], who argued that sheet structure and the shear stiffness between adjacent sheets will modulate changes in wall thickness. As argued in [46], this effect may be difficult to capture in traditional continuum mechanics models, because this phenomenon requires discontinuous strain fields. As for the Costa law parameters, variation in

*b*,

_{fs}*b*, etc., manifested clearly in their respective strain components (results not shown), but only the fibre stiffness parameter

_{fn}*b*affected many higher-level phenotypes.

_{ff}The overall (*a*) and fibre (*b _{ff}*) stiffness parameters dominated the variation in the whole-organ phenotype ‘change in volume’ (Δ

*V*). This result is in agreement with Stevens

*et al*. [22], who found that the maximal fibre strain parameter is the dominant factor for changes in volume in studies based on the pole-zero constitutive law [37]. On the other hand, variation in the endocardial fibre angle (

*α*

_{endo}) did not propagate to a change in volume (Δ

*V*) despite manifesting clearly in the long and short-axis diameter (Δ

*La*and Δ

*Sa*). The latter two effects tend to cancel out, as the long axis expands more when fibres are aligned more horizontally, whereas the opposite is true for the short axis.

Twist phenotypes stood out as being primarily governed by fibre angles (*α*_{endo} and *α*_{epi}) in a quite nonlinear fashion. Intermediate fibre angles tended to produce more twist than either flat or steep angles (simulation results not shown), possibly because oblique fibres contribute to torsion, whereas pure longitudinal and circumferential fibres do not. This nonlinear effect of microstructural parameters accounts for a lot of the variance that goes unexplained in the purely linear models (compare figure 4*b* with *c*).

### 3.3. Investigating interaction and nonlinearities in the genotype–phenotype map

A strength of the cGP approach is the ability to identify, e.g. dominance and epistatic effects by classical genetic methods such as variance decompositions, then drill down into underlying patterns of variability in the mathematical physiological model. Figures 5⇓–7 illustrate the interplay between microstructural and material parameters in determining spatial unevenness in stress and strain (the phenotypes with the largest residual variance in figure 4) for the different heart geometries. These are phenotypes of putative clinical relevance, partly because contraction can be stronger if fibre strain is high and evenly distributed at the onset of systole [23,47], but also because high stress and strain have been implicated in destructive remodelling [48].

The largest residual variance in figure 4*b* is observed for the spatial unevenness of fibre stress (var(*σ*_{ff})) in an eccentrically hypertrophic ventricle. As seen in figure 5, the mean effect of endocardial fibre angle on stress unevenness in this geometry is strongly non-monotonic, with stress becoming more unevenly distributed for both positive and negative deviations from the baseline of *α*_{endo} = 60°. For the other two geometries, there is little difference in unevenness between having one or two ‘alleles’ for steep *α*_{endo}. In genetic terms, the allele for steep *α*_{endo} shows *complete dominance* in the normal and concentric geometries, but *overdominance* in the eccentric case.

An uneven distribution of stress may be a driving force of microstructural remodelling in the heart [28,49]. For instance, figure 5 shows that the eccentric geometry results in a narrower range of fibre angles that minimize fibre stress. We hypothesize that ventricular hypertrophy may shift the optimum fibre angle and that the heart muscle fibres may reorient to even out the mechanical loads. This agrees with the observation of microstructural changes in dilated hearts [13].

Drilling down into what a one-dimensional summary of spatial unevenness showed above, figure 6 visualizes the raw spatial variation in fibre stress for the different heart geometries and fibre angles (keeping other genetic parameters at baseline). For *α*_{endo} = 30°, stress is highest in a thin layer near the endocardial surface, whereas steeper fibre angles distribute stress more radially. The one scenario that concentrates stress near the epicardial surface is the combination of steep fibre angle and eccentric geometry. This is in agreement with the finding in [23] that thinner heart walls tend to have highest stress further away from the inner surface.

Complementary to the above plots, figure 7 projects the simulated results onto two phenotypic dimensions (spatial mean and variability in fibre strain, *E _{ff}*) conditional on geometry and two ‘genetic’ parameters for material properties (fibre stiffness,

*b*) and microstructure (

_{ff}*α*

_{endo}). The clinical relevance is that elastic energy stored in diastole can contribute to an efficient systole if the fibre strain is high, but evenly distributed, using the Frank–Starling mechanism [27] of stronger contraction in response to higher strain. In our simulations, this is the case for the eccentric geometry at intermediate

*α*

_{endo}and low

*b*. Although eccentric hypertrophy primarily impairs contraction, which is not modelled in this study, we speculate that the high and even fibre strain may partially compensate for negative effects of eccentric hypertrophy. Figure 7 also illustrates differences in phenotypic variability among combinations of

_{ff}*α*

_{endo}and low

*b*, as well as a nonlinearity similar to that in figure 5. There was a positive phenotypic correlation between mean (

_{ff}*E*) and s.d. (standard deviation) (

_{ff}*E*) for parameter scenarios

_{ff}*within*the concentric geometry, although the correlation was negative

*between*geometries.

### 3.4. Towards causally cohesive genotype–phenotype modelling as a tool in personalized medicine

A major obstacle to achieving personalized medicine will be the ubiquitous context-dependence of genetic effects. Genetics offers many concepts for describing context-dependence at various levels, such as epistasis and genotype-by-environment interactions, but few tools for understanding the underlying mechanisms and including them in prediction models [5]. As demonstrated here and in earlier studies [39,50], cGP modelling is very well suited for integrating different types of variation, revealing context-dependent effects and identifying the underlying mechanisms.

A key step for application in personalized medicine will be new phenomics technologies and experimental programmes aimed at characterizing the genotype–parameter map in specific populations [5,51]. Furthermore, pragmatic application of cGP modelling must limit itself to a subset of the complex interplay of genetic, age-related, environmental and lifestyle factors. Effect modifiers outside the scope of the model proper may be included as statistical covariates or via phenotypic measurements that summarize the accumulated physiological effects, e.g. of lifestyle and environmental history (e.g. personalized heart geometries [52]). We believe that personalized medicine will benefit from a combination of cGP modelling and strategic phenotyping, equipping the digital patient with a combination of genetic sequencing information and context information such as familial history, occupational factors and phenotypic assays.

## Funding statement

This work was supported by the Research Council of Norway under the eVITA program, project number 178901/V30, and by the Virtual Physiological Rat Project, funded through NIH grant P50-GM094503. NOTUR, the Norwegian metacenter for computational science, provided computing resources under project nn4653k.

## Acknowledgements

We thank Joakim Sundnes, Simula, UiO for valuable comments on the manuscript.

- Received October 23, 2014.
- Accepted March 6, 2015.

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