Creating the Virtual Human

• Abstract
• Introduction
• Scales of the Modelling
• The Heart as an Example
  • Molecular level modelling
  • Cell level modellingm
  • Tissue level modelling
  • Whole organ modelling
  • Whole body modelling
• Discussion and Future Work Control
• Conclusion
• Acknowledgements
• References

Abstract
We report here on the progress made to date on developing an anatomically realistic, biophysically based model of the human body that incorporates the biochemistry, biophysics and anatomy of cells, tissues and whole organs. This work is the culmination of more than 25 years of research activity based at The University of Auckland. The comprehensive modelling framework being developed provides us with the opportunity to analyse integrative biological function in terms of underlying structure and molecular mechanisms.

Introduction
Modern medicine is currently benefiting both from the development of new genomic and proteomic techniques, based on our recently discovered knowledge of protein-encoding sequences in the human genome, and from the development of ever more sophisticated clinical imaging devices such as MRI, NMR, micro-CT, ultrasound imaging and optical tomography. This will mean that the clinical assessment of a patient’s medical condition could, in the near future, include information from both diagnostic imaging and DNA profile or protein expression data. To relate these two ends of the spectrum, however, will require very comprehensive integrative mathematical models of human physiology based on patient-specific quantitative descriptions of anatomical structures and models of biophysical processes which reach down to the genetic level.

We are thus on the verge of a model-based revolution in the biomedical field. Much publicity has been afforded to the completion of the first draft of the human genome sequence.^{[ 1 ]} However, rather than signal the end of an event, this achievement, remarkable as it is, really signifies the beginning of something much bigger, as the authors themselves allude to when they state, "In principle, the string of genetic bits holds long-sought secrets of human development, physiology and medicine. In practice, our ability to transform such information into understanding remains woefully inadequate". The challenge therefore is clear – to integrate this incredible wealth of information to allow the determination of structure and function at all levels of biological organisation.

This challenge has been recognised by the International Union of Physiological Sciences (IUPS), and confronted face on by the establishment of what has been termed the "Physiome Project".^[2,3] This concept was presented in a report from the Commission on Bioengineering in Physiology to the IUPS at the 32nd World Congress in Glasgow in 1993. The name comes from "physio-" (life) and "-ome" (as a whole), and is intended to provide a "quantitative description of physiological dynamics and functional behaviour of the intact organism". A satellite workshop, "On designing the Physiome Project", organised and chaired by the Chair of the IUPS Commission on Bioengineering in Physiology, was held in Petrodvoretz, Russia, following the 33rd World Congress in St Petersburg in 1997. A satellite meeting on the Physiome Project was held at the 34th World Congress of IUPS in Christchurch, New Zealand, in August 2001 and the Physiome Project was designated by the IUPS executive as a major focus for IUPS during the next decade.

The long-range goal of the Physiome Project is to understand and describe the human organism, its physiology and pathophysiology, and to use this understanding to improve human health. A major aim is to develop mathematical and computer models to integrate the observations from many laboratories into quantitative, self-consistent and comprehensive description. Integrated, biophysically based modelling is thus set for a revolution. Glimpses of this can already be seen. The pharmaceutical industry is beginning to appreciate the importance of mathematical modelling in reducing the massive cost of bringing a new drug to market and several companies have recently been formed to exploit biological modelling for drug discovery. The US Food and Drug Administration has indicated clearly that it would like submissions to include computer simulation and this in fact has been used by Roche to gain approval from the FDA for its cardiac calcium channel blocker Posicor.^{[ 4 ]}

In this paper we briefly describe the progress being made at The University of Auckland for handling the hierarchy of computational models and associated experimental data as we move towards the Physiome goal of the virtual human.

Scales of the Modelling
The wide range of spatial and temporal scales encompassed by the Physiome Project is shown in figure 1. It should be emphasised that no one model would encompass the 10⁹ dynamic range of spatial scales (from the 1nm pore size of an ion channel to the 1m scale of the human body) or 10¹⁵ dynamic range of temporal scales (from the 1 microsecond typical of Brownian motion to the 70 years or 109s typical of a human lifetime). Rather, it requires a hierarchy of models, such that the parameters of one model in the hierarchy can be understood in terms of the physics or chemistry of the model appropriate to the spatial or temporal scale at the level below. This hierarchy of models must range from gene networks, signal transduction pathways and stochastic models of single channels at the fine scale, up to systems of ordinary differential equations, representing cell level function, and partial differential equations, representing the continuum properties of tissues and organs, at the coarse scale. This key concept of being able to relate models across multiple spatial scales will be vital to the success of the Physiome Project, and is illustrated in figure 2.


Figure 1: Spatial (top) and temporal (bottom) scales encompassed by the Human Physiome Project.

The types of mathematical model appropriate to each spatial scale are also indicated. The last two images on the right in this figure are from anatomically based models developed by the Auckland Bioengineering group. The tissue image is a 3D confocal microscopy reconstruction of a transmural segment of rat heart by the group.

Figure 2: The process of integrating from cell (osteoclast) to tissue (trabecular bone) to organ (femur) to organ system (leg) is illustrated here.

The mechanical stress computed at the organ system level can then be fed back to the cellular processes controlling the balance of osteoblasts and osteoclasts in the bone-modelling unit.

The Heart as an Example
To demonstrate many of the modelling techniques and philosophies inherent in the Physiome concept, it is instructive to review an example. The methods demonstrated below have also been successfully applied to model the gastrointestinal, pulmonary and musculoskeletal systems and an image of our virtual human showing some of the organ model developments is given in figure 3. However, because of its relevance to mortality and morbidity, heart modelling is arguably the most advanced example of integrative biomedical modelling that currently exists, and we focus on that below (for more detail, see Hunter, Pullan & Smaill’s "Modelling total heart function"^{[ 5 ]}).

The heart is a variable speed biochemical pump. Electrical impulses originating in specialised cells and caused by the flow of ions across the membrane of the cell gives rise to a self perpetuating wave of activation that spreads through the entire myocardium. This activating wave initiates and, by the nature of its spreading, co-ordinates the contraction of the cardiac cells pumping blood to both the body and to the heart itself. This simple description belies the complexity of and the many sub-processes that exist within the beating of a heart at a large number of scales. The modelling approaches associated with each spatial scale are described below along with an outline of how each scale is integrated into a whole-heart-modelling framework which can be used to determine function in health and disease. 

Figure 3: The Auckland Bioengineering group’s virtual human, showing partial representations of the musculoskeletal, lung and gastrointestinal systems.

Molecular level modelling
At the molecular level protein chains interact to enable a cell to contract or regulate its internal or external ion concentration. Actin and myosin are long protein chains which progressively bind and unbind in a ratcheting motion to contract a cell along their length. Their function, as discussed above, was first elucidated by Andrew Huxley^{[ 6 ]} using a mathematical model of two, coupled, partial differential equations which account for the spatial strain and the rates of attachment and detachment. This approach has since been extended in a number of models to include a larger number of binding states between the two proteins and a related increase in the number of coupled equations.

The interaction of proteins also controls the flux of ions into or out of a cell by lining pores in the cellular membrane. There are a large number of models which use a single ordinary differential equations or a system of them to quantify how the overall conductance of the channel changes over time or is affected by the concentration of ions which may bind to sites on the channel.

The structure and function of the pore proteins are determined by their encoded gene sequence and their structure. A mutated gene can affect ion motility across the cell membrane, affecting the activation and recovery of the cell, which, in the worse case scenario, can lead to cardiac arrhythmia and death. The molecular structure of potassium (K) channels, while still controversial, is likely to be finalised shortly from xray crystalography, and the finalisation of sodium (Na) and calcium (Ca) channels will follow. The determination of function using molecular modelling based on the xray and NMR structural data is still extremely difficult but has at least had success with open K channel conductance. The voltage dependence is still being determined and gating kinetics are beyond current computational reach at the moment, but higher performance computers and better coarse-grained model approximations will bridge the gap from molecular structure to Hodgkin-Huxley kinetic models within the next few years.

Cell level modelling
Cellular models which describe the chemical reactions within the cell are used to determine functions such as excitability, tension generation and metabolic energy production. Each of these models is expressed as a series of coupled rate equations where parameters include temperature, dimensions and ion concentrations of the cell.

Well-known examples of such models include Noble 98^{[ 7 ]} and LRII^{[ 8 ]} for the electrical activity, HMT^{[ 9 ]} for mechanics and the Salem model^[10] of cardiac metabolism. Each of these models can be defined in three parts:

1. model parameters, such as physical constants and fixed reaction rates or channel conductances;
2. state variables, which define the state of the cellular system at any given point in time such as ion concentrations or channel gating states; and
3. derived variables, such as ion current through the membrane.

Typically, functions of the derived variables and model parameters are used to construct expressions which determine the rate of change of each state variable over time. This system of equations is then integrated through time to determine transient cellular function. One difficulty of coupling distinct cellular functions together is the very different time scales of each system. For example, at the cellular level excitation occurs on the scale of milliseconds whereas metabolic processes can take many minutes to develop fully. This problem is often overcome using mathematical techniques which approximate a fast system with a fine time scale and a slow system with a coarse time scale while still storing the values of state variable at all times. Recent coupled cellular models encompassing excitation, contraction and metabolism contain up to 50 state variables, 100 parameters and 80 derived variables.

Tissue level modelling
Cardiac structure plays a pivotal role in the contraction of the heart – cells individually contract along their length but this cell shortening in itself is not sufficient to generate the volume changes seen in a whole heart beat. Significant "sliding" of groups of cells over other groups also occurs, the directions of which are determined by the arrangement of the cells within the heart, and the distribution of collagen. Conduction speeds through tissue are also critically dependent on tissue structure being, typically, ten times faster in the direction orthogonal to the sheets.

To model both health and disease these effects of tissue structure are clearly beyond the scope of single cell models. Thus, to determine tissue and, ultimately, whole organ function, a modelling framework to which structural information can be added to cellular function is required. The cell models (outlined above) are used to provide sources of current, tension and oxygen tension. Using a finite element mesh to spatially represent the tissue sample, a cell model is embedded at evenly spaced grid points where each grid point is treated as a "black box’’, the source/sink characteristics of which are determined by the complex ion kinetics of the underlying cellular model.

The spread of current is then modelled by numerically solving the well-known bidomain equations underlying myocardial activation and extracellular current flow which are formulated with a diffusion tensor based on the structural model. Within the tissue model, conductivity and capacitance determine the electrical coupling between approximation points and ultimately the spread of activation throughout the tissue. The anisotropic nature of cardiac tissue means this conductive coupling is represented by a tensor with values that are determined from the tissue micro-structure that has its local axes aligned with the fibre, sheet and sheet-normal directions in the tissue.

In addition to the active tension generated by the embedded cellular models, the mechanical behaviour of tissue models is dependent on both the passive stress strain relationships. This passive relationship is quantified using non-linear constitutive law which determines the tissue stiffness locally about the three axes defined by the tissue organisation (fibre, sheet and sheet-normal). The parameters of this constitutive law are estimated from biaxial tension, compression and shear tests on small samples of ventricular tissue. Using these material properties, tissue deformation is determined using the finite element method to solve the governing equations of finite deformation elasticity.

Figure 4: This figure illustrates schematically some of the major processes involved in the electrical, contraction and coupled electro-mechanics simulation on a two dimensional tissue block demonstrating the propagation of a spiral wave.

Whole organ modelling
The step up in scale from tissue to whole organ modelling is initially one of adding geometry and characterising the regional variation in structure. Careful dissections of hearts^[11] together with detailed imaging^[12] have gathered significant information on the unique organisation of cardiac myocytes within the heart. These detailed measurements have been encapsulated within whole organ continuum models using non-linear optimisation fitting techniques and piecewise parametric modelling. Figure 5(a) shows a finite element mesh constructed from fitting data points from pig ventricles.

Figure 5: Figure 5(a) shows a finite element mesh fitted from measurement of the cardiac ventricles. Figure 5(b) shows the regional variation in principle stresses and their direction at end-diastole at midwall points of the finite element mesh.

Using the tissue modelling techniques applied to this anatomical model, the dynamics of heart excitation and contraction can be determined at the whole organ level. Figure 5(b) shows the directions and magnitudes of the principle stresses in the heart at the end of contraction, illustrating the ability of a model to calculate and display data which is extremely difficult to measure through experiments during the heart cycle.

To fully understand heart function in health and disease, a model of coronary blood flow needs to be added which can predict the regional distribution of blood flow through the heart wall. The anatomy of the coronary vasculature is measured on the epicardium and fitted to the surface of the heart model. To represent smaller vessels which branch into the heart wall a nonlinear network growth algorithm is used to add spatial information to detailed topological measurements made from silicon casts of the coronary vessel.^[13] Figure 6 shows the resulting anatomically based finite element model embedded in ventricular geometry.^[14,15]

A prediction of the regional and temporal variation in coronary blood flow is achieved by applying computational fluid dynamics methods on the domain of the anatomically generated model of the large coronary vessels.

Figure 6: An anatomically based model of the largest six generations of the coronary arterial network of blood vessels embedded in cardiac geometry.

Whole body modelling
Physically, cardiac electrical activity is connected to current flow in the torso via the extracellular potential field. Cutaneous electrodes are routinely used in ECG clinics all around the world to record such activity. Mathematical models of human torsos have been developed which can incorporate the information from the cell and organ level models described above, and thus simulate ECG signals. While the individual details of a single cell’s electrical activity cannot be uniquely distinguished at the body surface, to accurately simulate an ECG still requires much of the detail of the cellular level models described above. An example of the torso potentials generated from one such model is illustrated in figure 7.

Figure 7: The body surface potentials calculated at 25 ms using an anatomically accurate torso model.

Voltage is scaled between -0.1 mV and 0.2 mV and each band represents a 0.025 mV increment.

Discussion and Future Work
The next step in this project is to extend the results outlined above to mechanisms of disease and how they relate to heart function. This progression relies on the geometric models of the heart and torso presented above to be combined with detailed cellular and, potentially, genetically controlled sub-cellular models. This would result in a considerable increase in computational time. Additional organ scale coupling will need to be added between the predicted distribution of coronary blood flow and our model of cardiac metabolism such that the supply of oxygen determines the nature of energy production and the resulting effects on myocyte production. Ultimately, we envision an accurate beating heart in which patient-specific profiling (using genetic as well as anatomical data) may be incorporated, opening the way for personalised medicine prescribed in response to model-based predictions.

Conclusion
We have outlined a biophysically and anatomically based computational framework for simulated integrated whole organ cardiac function which provides the potential to understand the manifestations of heart function at the cell, whole organ and torso spatial scales. Such integrative biophysically based modelling is a central theme of the worldwide Physiome Project, which aims to use modelling to understand and describe the human organism and ultimately use this understanding to improve human health.

Acknowledgements
The contributions to the work described in this paper by the many past and present members of the Bioengineering Institute are warmly acknowledged. The work has also benefited from grants from a number of funding agencies including the Royal Society of New Zealand Marsden Fund and the Health Research Council of New Zealand.

References

1. International Human Genome Mapping Consortium. A physical map of the human genome, Nature 2001;409:934–941.
2. Hunter PJ, Borg TK. Integration from proteins to organs: the Physiome Project. Nat Rev Mol Cell Biol 2003; 4:237–243.
3. Crampin EJ, Halstead M, Hunter PJ, Nielsen PMF, Noble D, et al. Computational physiology and the physiome project. Exp Physiol 2004 89:1–26.
4. Noble D, Colatsky TJ. A return to rational drug discovery: computer-based models of cells, organs and systems in drug target identification. Emerging Therapeutic Targets 2000;4:39–49.
5. Hunter PJ, Pullan AJ, Smaill BH. Modelling total heart function. Annu Rev Biomed Eng 2003:5: 147–177.
6. Huxley A F. Muscle structure and theories of contraction. Prog Biophys Biophys Chem 1957;7:255–318.
7. Noble D, Varghese A, Kohl P, Noble P. Improved guinea pig ventricular cell model incorporating a diadic space, IKr and IKs, and length- and tension-dependent processes, Can J Cardiol 1998 Jan;14(1):123–34.
8. Luo C, Rudy Y. A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes. Circ Res 1994 Jun;74(6):1071–96.
9. Hunter PJ, McCulloch AD, ter Keurs HE. Modelling the mechanical properties of cardiac muscle. Prog Biophys Mol Biol 1998;69(2–3):289–331.
10. Salem JE, Saidel GM, Stanley WC, Cabrera ME. Mechanistic model of myocardial energy metabolism under normal and ischemic conditions. Ann Biomed Eng 2002 Feb;30(2):202–16.
11. LeGrice IJ, Hunter PJ, Smaill BH. Laminar structure of the heart: a mathematical model. Am J Physiol 1997;272(5 Pt 2):H2466–76.
12. Hooks DA, Tomlinson KA, Marsden SG, LeGrice IJ, Smaill BH, et al. Cardiac microstructure: implications for electrical propagation and defibrillation in the heart. Circ Res 2002 Aug 23;91(4):331–8.
13. Kassab GS, Rider CA, Tang NJ, Fung YC. Morphometry of pig coronary arterial trees. Am J Physiol 1993 Jul;265(1 Pt 2):H350–65.
14. Smith NP, Pullan AJ, Hunter PJ. Generation of an anatomically based geometric coronary model. Ann Biomed Eng 2000;28(1):14–25.
15. Smith NP, Pullan AJ, Hunter PJ. An anatomically based model of coronary blood flow and myocardial mechanics. SIAM J Appl Math 2002;62:990–1018.