Multi-scale computational modelling in biology and physiology

Abstract

Recent advances in biotechnology and the availability of ever more powerful computers have led to the formulation of increasingly complex models at all levels of biology. One of the main aims of systems biology is to couple these together to produce integrated models across multiple spatial scales and physical processes. In this review, we formulate a definition of multi-scale in terms of levels of biological organisation and describe the types of model that are found at each level. Key issues that arise in trying to formulate and solve multi-scale and multi-physics models are considered and examples of how these issues have been addressed are given for two of the more mature fields in computational biology: the molecular dynamics of ion channels and cardiac modelling. As even more complex models are developed over the coming few years, it will be necessary to develop new methods to model them (in particular in coupling across the interface between stochastic and deterministic processes) and new techniques will be required to compute their solutions efficiently on massively parallel computers. We outline how we envisage these developments occurring.

Fig. 1

Levels of biological organisation.

Fig. 2

Cutaway diagram of a typical MD potassium channel set-up. Two subunits of the tetramer have been removed for clarity. Shown are the lipid bilayer in a box of water, two subunits of the pore region of the protein, two ions in the filter separated by a water molecule, some water in the cavity and the waters positioned just outside the filter region.

Fig. 3

Effect of regional differences in action potential shape and duration on repolarisation and the electrocardiogram in the canine heart. Top panel left, action potentials for endocardial, mid-myocardial and epicardial cells during pacing at 1000 ms intervals. Top panel right, APD restitution for the three cell types obtained in a thin strip of simulated tissue using an S1 S2 protocol, with a basic (S1) cycle length of 500 ms. The cell model used is the four-variable reduced model described in Cherry and Fenton (2004) with parameter sets fitted to experimental recordings from canine endocardial, mid-myocardial and epicardial cells as described in Cherry et al. (2003). The parameter sets were obtained by personal communication from E.M. Cherry and F.H. Fenton. Lower panel shows colour-coded repolarisation times in long and short axis slices through the Auckland canine ventricle geometry (Nielsen et al., 1991), obtained using an explicit FD solution of the monodomain model with a space step of 0.25 mm, a time step of 0.05 ms, and anisotropic diffusion with coefficients of 0.002 cm² ms^–1 and 0.0005 cm² ms^–1 along fibres and transverse to fibres, respectively. The four images show repolarisation patterns during steady pacing at 500 ms intervals from an epicardial distribution of sites of earliest activation based on experimental data (Scher and Spach, 1979). The top two images show repolarisation with uniform epicardial cells, and the lower two images show repolarisation patterns obtained during steady pacing with layers of endocardial, mid-myocardial and epicardial cells, each of equal thickness. The right-hand panel shows a simulated electrocardiogram assuming the heart is immersed in an infinite volume conductor. The traces shown were computed for the uniform epicardial model (blue) and for the heterogeneous model (red), showing that the different pattern of repolarisation associated with heterogeneity produces an upright instead of a biphasic T wave. Courtesy of Dr. Richard Clayton.

Fig. 4

A rabbit ventricular geometry model (Vetter and McCulloch, 1998) with smooth epicardial and endocardial surfaces and anatomically realistic fibre orientation was discretised using an unstructured grid with an average discretisation of 250 μm. The bidomain equations were discretised on this non-uniform grid using linear tetrahedral elements with a time step of 8 μs. The active membrane behaviour was described by the rabbit ventricular Puglisi model (Puglisi and Bers, 2001) incorporating an electroporation current (De Bruin and Krassowka, 1998) and a hypothetical I_a current (Cheng et al., 1999). Details of the applied numerical methods are found in Plank et al. (2007). Two plate electrodes, a stimulation electrode and a grounding electrode, were used to stimulate the ventricles by delivering a train of ten pulses (see top panel). Subsequently, electric activity was simulated for another 2 s. Simulations were carried out with C_m=1 μF cm^–2 and β=1400 cm^–1. Initially, conductivity along the fibres was set to σ_il=1.74 mS cm^–1 and σ_el=6.25 mS cm^–1, and transverse to the fibres to σ_it=0.19 mS cm^–1 and σ_et=2.36 mS cm^–1, in the intracellular and interstitial domain, respectively (Clerc, 1976). The conductivity of the surrounding fluid was set to σ_b=1.0 mS cm^–1. To determine suitable parameters which lead to a re-entry under the given protocol, the basic cycle length (BCL) of the pulses and the wave length of the tissue, λ, were varied. From the series of simulations, a standard was chosen with a BCL of 200 ms and λ reduced to 0.66 of the nominal λ, as computed with the default conductivity settings, leading to a sustained figure-of-eight re-entry circulating around the apex (see bottom panels). Courtesy of Dr. Gernot Plank.