Computational multiscale modeling in the IUPS Physiome Project: Modeling cardiac electromechanics

Mathematical modeling may be used to integrate data from electrical and mechanical cardiac experiments in order to test hypotheses that are concerned with multiple spatiotemporal scales and functions. This kind of integration across a variety of size and time scales is among the most advanced examples of organ system modeling [1–3]. A goal of the IUPS (International Union of Physiological Sciences, www.iups.org) Physiome Project (www.physiome.org.nz) is to develop the technology and methods required to simulate the behavior of biological organisms via numerical simulations using computational and mathematical models. Such simulations require the integration of multiple types of physics over a wide variety of spatial and temporal ranges, for example, spatial scales of 10⁻⁹ meters for subcellular structures up to approximately one meter for the human body, and molecular events occurring on the 10⁻⁶-second time scale up to the human lifetime of the order of 10⁹ seconds (Figure 1).

Figure 1

As the complexity of computational models increases, formal vocabularies are needed to reduce the growing heterogeneity of biological and mathematical expressions. Standards are being developed to formalize the description of experimental data and mathematical models of physiological processes [4, 5]. Ontologies that incorporate semantic descriptions of modeling concepts eliminate ambiguities in the modeling environment. To this end, ontologies and representation languages that include these ontologies are being developed under the IUPS Physiome Project, which facilitates communication of biological models to researchers through tools for building, sharing, interpreting, and visualizing models [6].

With the development of these languages comes the ability to create and populate repositories of models that are freely available for use by the scientific community. Currently, the most advanced of the “physiome” representation languages is CellML, a markup language based on open standards. Initially designed for application to models of cellular electrophysiology and reaction pathway models, CellML has since been used in a wide range of mathematical models, including constitutive laws for continuum mechanics. The term physiome refers to the quantitative and integrated description of the functional behavior of the physiological state of an individual or species.

In this paper, we review the development of a computational modeling framework that enables scientists to perform numerical simulations using models of integrative physiological function across the cellular, tissue, and organ spatial scales. The framework was initially developed with a focus on overcoming the problem of modeling the tightly coupled electromechanical function of the heart, but the implementation is sufficiently general that applications to other areas of modeling are underway.

The framework we have developed allows the use of CellML for the specification of model- and simulation-specific mathematical equations—for example, cellular electrophysiology models and passive mechanical functions. The simulation framework is built upon the extensive base provided by the CMISS software package developed at the Bioengineering Institute at The University of Auckland, New Zealand. CMISS is an interactive computer program for Continuum Mechanics, Image analysis, Signal processing and System identification. It is also a modeling environment that allows the application of several mathematical techniques to a variety of complex bioengineering problems.

CellML (www.cellml.org) is an XML-based language developed by the Bioengineering Institute at The University of Auckland [4] (originally in collaboration with Physiome Sciences Inc., but now entirely supported by New Zealand's Public Good Science Funding). CellML is a language designed to store and exchange computer-based biological models; where appropriate, the language builds upon existing XML standards such as MathML (www.w3.org/Math) for the specification of mathematical equations and the Resource Description Framework (www.w3.org/RDF) for the encapsulation of metadata.

CellML provides a relatively basic set of tags that can be used to mark up complex interactions between a set of mathematical equations represented in the MathML language. See Figure 2 for an example of CellML code. Although CellML clearly provides a means for making models available to researchers for validation and study, these models must also be published and peer-reviewed before being accepted by the modeling community [7].

Figure 2

Through the use of CellML and the open-source tools that are now becoming available (cellml.sourceforge.net), an author of a model is able to describe a model in such a way that others are able to incorporate it into their own computational codes or simulation package of choice. The model equations can be specified just once and used in all implementations and publications of the model. Similarly, all boundary and initial conditions required for a particular computational experiment can be specified just once and used by the community. With the establishment of a publicly accessible repository of models and simulation tools, authors are able to submit validated models and simulation results to ensure that other investigators are able to accurately reproduce their simulations. In this paradigm, the onus of model validity no longer rests with the model user, but with the model author and the software engineers who implement the CellML application libraries and tools. The repository of models and simulation results consists of data that the author asserts is an accurate representation of a model; it provides a much larger basis for testing the code and ensuring robustness and compatibility.

As model repositories are being developed, researchers need to be confident that a given representation or implementation of a model is “accurate”—that is, the degree to which the representation or implementation reflects the reference description of the model, or the accuracy with which the underlying phenomena represented by the model are reproduced. Standards are being developed that will address these concerns [8] and will be incorporated into CellML models through the use of curation metadata that will be crucial for the general acceptance and use of a CellML model repository.

The CellML model repository contains models from peer-reviewed journal publications. When researchers make a CellML version of a published model available in the repository, this does not guarantee that the model is error-free, because the original publication may have errors, and the corresponding CellML model versions, for issues of provenance, must faithfully reproduce the mathematics in the paper, errors included. However, it is clearly desirable to create a second version of the CellML-encoded model in which various checks have been carried out, including 1) checking that the defined units of measurement are consistent; 2) checking that all parameters and initial conditions are defined; and ultimately 3) checking that running the model reproduces the results published in the paper. As the CellML standard and the model repository gain acceptance by journals, it should be possible to work with authors to achieve this refined version of the model at the time of publication. For example, such refinement has happened during publication of a metabolic model [9]. A further level of model curation is anticipated in which models are checked for the extent to which they satisfy physical constraints such as conservation of mass, momentum, and charge in chemical and physical interactions.

Owing to the completely generic specification of the language, CellML has a much broader range of applicability than the name suggests. CellML is considered to be a language suitable for the description of annotated mathematics, i.e., the specification of mathematical equations using MathML, the interrelationships between equations, and the connections between variables contained in the equations. Model authors and users are further able to annotate models by using metadata associated with any data contained in the model.

Other XML-based languages have been developed to aid the exchange of mathematical models within various scientific communities. These languages have tended to be tied to specific domains in terms of their use and the actual definition of the language syntax. For example, CellML is often compared to the Systems Biology Markup Language (SBML), which is currently specifically for use with biochemical network models, although future development plans for the language begin to introduce concepts similar to those in CellML. A particularly powerful feature introduced in CellML 1.1 is the ability to reuse models, or parts of models, by “importing” the relevant mathematical equations and variable definitions into a new model. As we have suggested, the generic syntax used in the CellML language implies no specific domain of application, which allows its use in a wide variety of scenarios that allow model authors to annotate models with domain-specific metadata.

As mentioned, CMISS (www.cmiss.org) is a computational package for modeling the structure and function of biological systems. In particular, it is designed to model the anatomy and behavior of organ systems (e.g., cardiovascular, respiratory, and special sense organs) from the component organs (e.g., heart, lungs, and eyes), while also considering the cellular and subcellular scales and the coupling that occurs between and within all of these levels. Equations derived from physical laws of conservation, such as conservation of mass, momentum, and charge, are solved in order to predict the integrative behavior of an organ, given descriptions of the anatomical structure and tissue properties. The tissue properties used in these organ simulations can incorporate tissue structure and cellular processes, together with spatial variation of the parameters, such as parameters associated with mechanical compliance and electrical conductivity, that characterize these processes. CMISS has facilities for fitting models to geometric data derived from various imaging modalities (e.g., MRI, CT, and ultrasound) and has a rich set of tools that permit graphical interaction with the models and display of simulation results. A graphical user interface, an interactive console-based interpreter, or batch-mode scripting may be used to control CMISS.

CMISS originated in the doctoral work of Peter Hunter [10] as a finite-element program for stress analysis of large deformations in the heart. The package has since evolved into a general-purpose biological-systems modeling tool, used in the areas of continuum mechanics, image analysis, signal processing, and system identification. Recently, work has begun to modularize CMISS in order to enable the development of specialized and focused tools for various medical and other applications. The main academic goal of CMISS is to support the IUPS Physiome Project. See [6] for a detailed review of the abilities of CMISS in relation to the heart.

As the Physiome Project evolves, it is essential to provide both programmer and user access to the technologies being developed. For example, application developers must be able to access the model repositories and the data contained within model-representation documents, while users must be able to interact with specific models and perform simulations. To achieve this support, the user interface of CMISS is separated from the main application and is released under an open-source license. As part of this software evolution, standard application program interfaces (APIs) are being developed for various components of the software, which is being divided into separate modules that can be linked into external applications.

With the development of these modules, customized user interfaces can be created for specific modeling or simulation platforms along with interfaces capable of browsing model repositories. To aid in the sharing of information, most of these interfaces will be “Web-deliverable” either through the use of standard web technology or via a Mozilla extension that implements the APIs of the open-source CMISS modules. Mozilla extensions are applications that can be added into an existing Mozilla-based web browser (e.g., Firefox** or Mozilla) to provide extended functionality. Using this extension, graphical interfaces to Physiome technology can be specified using the XML user interface language (XUL) and delivered via the Internet to any XUL-enabled client such as the Firefox web browser.

Our goal was to implement a process for enabling the specification of model-specific mathematical equations in CMISS through CellML. For example, prior to this work, a new cellular model in CMISS was implemented by writing the equations in Fortran and adding this to the CMISS code base. This method has two main disadvantages: The author of the Fortran code is usually translating a published piece of work, which can lead to human errors in the translation, and the generated code is very specific to CMISS and thus not easily transferred to another modeling or simulation package. This also assumes that the published model itself is free from errors. Like problems associated with cellular models, the hard-coding of other mathematical equations into the CMISS code base has disadvantages.

To avoid these problems, CMISS is capable of importing mathematical models from CellML (Figure 3). This provides the ability to store and simulate models in an open standard, even though the models may conceivably originate from various sources.

Figure 3

The first step in implementing this ability in CMISS required us to develop a standard API for use when accessing a CellML model description. The most recent API implementations are freely available from the CellML website, cellml.sourceforge.net. With an API defined, the ability to import CellML was added to CMISS, allowing the definition of mathematical models via CellML. This also required the capability of translating the mathematical expressions from the CellML model (stored as MathML) into a dynamically loadable object that can be utilized by CMISS during a given simulation. Given the structure defined by MathML, the translation to such an object is reasonably straightforward (Figure 4).

Figure 4

Using the described approach, a scientist can develop a cellular model using software that is well suited to single-cell modeling or that may be used in experimental work. If the software is capable of exporting CellML, it is possible to use a cellular model and perform tissue and whole heart simulations that are based on the model.

In order to import CellML models into tissue representations, we must be able to specify spatial variations of these models and their parameters within the larger-scale model. For example, when modeling the spread of excitation from the pacemaker cells in the sino-atrial node into the atria, a modeler would typically use different cellular electrophysiological models for each of these two regions of tissue. The variation of channel distributions through the ventricular wall of the heart (Figure 5) is another example of the need for spatial variation. Here, the modeler may want the same equations represented at all points in the ventricles, but needs to specify spatial distributions of channel densities.

Figure 5

The modeling framework that we have developed uses cellular models of cardiac electromechanics to drive the dynamic functional behavior of tissue and organ continuum models. Electrical excitation and mechanical deformation at the tissue and organ spatial scales, in turn, modulate changes in cellular model material properties.

Cardiac cellular electrophysiological modeling is a well-established field with numerous existing models that cover a large range of species and cellular phenotypes (see [7, 11, 12] for recent reviews). Similarly, the mechanical behavior of cardiac cells has been extensively investigated through the use of mathematical models. In recent years, modeling of cellular metabolism and energetics [13, 14] has been considered an established field, and this modeling has gained prominence because of its relevance to many dysfunctions of the heart. Each of these functions is closely integrated within a cardiac cell, and some models have been developed that reflect this tight coupling at the cellular or subcellular level [13, 15]. Another approach requires the development of methods that enable independent models to be coupled together, and this approach has been shown to work particularly well when coupling cellular electrophysiological and mechanical models [16].

Cellular models have been developed with varying levels of biophysical detail. For example, when modeling various phenomena, different approximations can be made to simplify the cellular model in order to dramatically decrease the computational demands required for simulations using these models. Thus, cellular models are typically identified as either “low-dimensional” or “biophysically based” models. Low-dimensional models are used to represent gross behavior of the system, sometimes based on the use of mathematical analysis to reduce complex models, given certain constraints. A classical low-dimensional model of the cellular action potential is the FitzHugh–Nagumo model [17, 18], based on a phase-plane analysis of the Hodgkin–Huxley nerve axon model [19] and modified in several cardiac-specific action potential models [20, 21].

In contrast to the relatively simple equations of low-dimensional models, biophysical cellular models attempt to accurately represent detailed physiological processes and mechanisms that underlie the phenomena being modeled. In the case of cardiac electrophysiology at the whole-cell spatial scale, this includes the dynamics of various ionic species and the gating kinetics of various proteins to permit or block the transport of ions between distinct compartments. Such models generally consist of large systems of stiff ordinary differential equations. Models based on the work of C. Luo and Y. Rudy [22, 23] and D. Noble [24] have been widely used and adapted to various specific situations. As the quantity and quality of experimental data and techniques improve, greater biophysical detail can be extracted. Furthermore, models are now being developed to reproduce the effects of genetic mutations that govern the dynamics of specific transmembrane ion channels. These models are beginning to incorporate more realistic stochastic behavior of protein populations contained in single cells or populations of entire cells, leading to significant increases in the complexity of the models [15, 25, 26].

Modeling cardiac electromechanics on a spatial scale larger than for the single-cell models described above requires the coupling of two processes: the spread of electrical excitation through the tissue and the mechanical response of the tissue. In a normal mammalian heart, the spreading wave of electrical excitation triggers contraction of the cardiac muscle, which is responsible for the pumping blood. While the aforementioned two processes can be modeled independently, it is well established that mechanisms of excitation–contraction coupling and mechano-electrical feedback are tightly linked [27].

A large body of research exists for modeling the spread of electrical excitation throughout cardiac tissue as well as the whole heart [12, 28]. Continuum models are based on the assumption that the length scales of the physically observable phenomena are large in comparison to the underlying discrete structure of the material. Our group has developed numerical simulation techniques for solving continuum models of electrical excitation, based either on the bidomain model (i.e., a model involving solving for intra- and extra-cellular potential fields, assuming two interpenetrating domains) [29, 30] or an eikonal-type model for activation times, which involves instantaneous solving for electrical activation times throughout the solution domain [31].

Although we have the tools to model the full bidomain model, in this work we have reduced the complexity of the simulations by using the simplified monodomain model [32] and neglecting the effect of extracellular potential on the electromechanical behavior of the tissue. These assumptions are valid when simulating the normal spread of electrical excitation, but the full bidomain model would be required in other circumstances, such as when defibrillation shocks are being applied [33]. Complexity reduction allows a significant saving in terms of computational cost in both memory requirements and solution times.

Similarly, a large body of research also exists for modeling the mechanical behavior of cardiac tissue (see [28] and [34] for reviews). In the past, many models of cardiac mechanics considered only the passive properties of the muscle for reasons of model complexity and the lack of experimental data during the systolic phase of the cardiac cycle (i.e., during the contraction of isolated tissue preparations). This allowed quasi-static models of finite-deformation elasticity to be applied to the heart with great success. Finite-element continuum models of cardiac mechanics are the most prevalent in this field, and the use of high-order interpolation of fields, in which our group specializes, is well suited to these types of models.

The development of models of electrophysiology and mechanics has largely occurred independently, and only recently have we begun to obtain the computational power and experimental data required to develop models of electromechanics in cardiac tissue. Various approaches have been used in the development of cardiac electromechanics models providing varying levels of physiological detail and interactions between electrical and mechanical processes. Such tissue models consider tight interaction between mechanics and electrophysiology using low-dimensional cellular models [35–37], and also models with less interaction between the mechanics and electrophysiology during a simulation (e.g., either excitation–contraction coupling or mechano-electrical feedback). These latter models are based on more biophysically detailed cellular models (e.g., [38–40]). In this case, models have typically solved for electrical activation times and used these times to trigger local active contraction of the cardiac tissue. The activation times can be computed either from a simulation of electrical activation or through the use of an eikonal model to solve for activation times directly [31]. In these loosely coupled frameworks, the spread of electrical activation is calculated independently of both mechanical deformation and mechano-electrical feedback mechanisms such as stretch-activated channels and calcium buffering by contractile proteins.

In the work described in the current paper, we have used a large-scale, high-order-interpolation, finite-element-based method for solving mechanics, coupled to a small-scale, low-order interpolation method for solving electrical activation in order to produce a technique for the numerical solution of biophysically detailed cardiac electromechanics models [41]. Our tightly coupled framework uses cellular models of electromechanics to drive the dynamic functional behavior of the model, while the properties of the cellular models are modulated by both the electrical excitation and the deforming mechanical model.

Owing to the computational resources required for modeling three-dimensional electromechanics, we have been limited to simulations that use either complex cellular models with simplified (small) geometrical models or low-dimensional cellular models with more anatomically based (large) geometries. Here we present results obtained from a cube of tissue, and some preliminary results from a simplified geometric model of the cardiac left ventricle, as an illustration of the application of the framework discussed above. Further analysis and discussion of the left-ventricular results can be found elsewhere [41].

We first present simulation results from the cube shown in Figure 6. We performed simulations to investigate the response of this tissue block to the applied electrical stimulus using two cellular electromechanics models: a low-dimensional and a biophysical model. The low-dimensional model was obtained through the coupling of the Fenton–Karma (FK) cardiac action potential model [42] to the Hunter–McCulloch–ter Keurs (HMT) mechanics model [43], following the procedure described in [16]. For the biophysical model, we coupled the HMT mechanics model to the most recent of the Luo–Rudy-based models developed to investigate various genetic mutations [25, 26, 44], with the addition of a more biophysically detailed model of calcium dynamics [45]. In the following, we denote the low-dimensional model as the FK–HMT model and the biophysical model as the N–LRd–HMT model.

Figure 6

Figure 7 illustrates a summary of the results from the biophysical simulation using the N–LRd–HMT cellular model for the tissue cube shown in Figure 6. In Figures 7(a) and 7(b), the wave of electrical activation can be seen advancing through the cube from the stimulus (isoelectric surfaces are shown in color, with blue and red shading respectively representing the most negative and the most positive potentials). Following activation, the adjacent tissue contracts along the fiber direction [shown in Figure 6(a)] due to the dynamic development of tension within the cells triggered by the electrical excitation (the axis of cardiac cells is aligned with the fiber direction). Given the incompressible nature of cardiac tissue, expansion occurs in the other two dimensions. Following excitation, the tissue recovers and returns to the initial resting state.

Figure 7

Figure 7(c) shows key transients from the cellular models at the spatial locations within the tissue cube shown in Figure 6(b). The transients shown include the cellular action potentials, dynamic tension, and a measure of cellular length. These transients show the smooth propagation of electrical excitation through the tissue, from the cyan point on the stimulated face to the red point farthest from the stimulus. Following a time delay (where electrical excitation triggers internal cellular processes that give rise to tension generation), the dynamic tensions of the cells begin to rise, causing the cells to shorten and hence causing the tissue to contract. Also illustrated is the initial stretch of tissue in regions distal to the stimulus due to the passive cells being stretched as the actively contracting tissue acts against the applied displacement boundary conditions, as described in Figure 6. As tissue is activated, it follows a similar contraction transient.

The dip into negative tension shown in Figure 7(c) was unexpected and was initially thought to result from poor numerical convergence. Further investigation revealed that this was not the case, and that the negative dip could be attributed to mismatched material parameters and deficiencies in the model. Material parameters for the cellular, electrical activation, and finite elasticity components of the model have been taken from previous studies in which the parameter values of each component have largely been established independently. Further material parameter estimation studies are required to determine a more appropriate set of parameters in order to better match experimental observations. Such studies are currently expensive to perform because of the computational requirements of these simulations (Figure 8). In addition to material parameter mismatches, acknowledged deficiencies exist in the active contraction model that contribute to non-physiological behavior in both the initial phase of the tension transient and during relaxation. These deficiencies have been addressed at the cellular level in a recent study [46].

Figure 8

Both the FK–HMT and N–LRd–HMT cube models consist of 64 high-order finite elements used for the solution of the equations of finite deformation elasticity and 46,656 low-order finite elements for the electrical propagation solution representing 35,937 cells at a spacing of 0.125 mm. The simulations were performed on an IBM eServer* p690 computer (frequently referred to as “Regatta”) with 32 1.3-GHz POWER4* processors and 32 GB of main memory.

The simulation using the FK–HMT cellular model required approximately 550 MB of memory while the N–LRd–HMT model required 1.5 GB, which will be of concern for model extension to more realistic geometries that require millions of cells. CPU and wall clock times for these simulations are summarized in Figure 8, with the simulations performed using eight processors and the shared-memory parallel implementation of CMISS. In Figure 8, the simulation time for each of the models is split into the three main components: solution of the electrical activation model, solution of the finite elasticity mechanics model, and the update step transferring data between the two models to integrate the feedback effects.

As shown in Figure 8, the main difference in solution time between the FK–HMT and N–LRd–HMT models is due to the electrical propagation model. In this step, the full cellular model is integrated over a time step at each of the 35,937 cells in the tissue cube, and the extra complexity in the biophysical N–LRd–HMT cellular model significantly increases the computational time required for model solution. During the mechanics solution step, no temporal integration is required in order to update the steady-state active tension solution, so we simply evaluate the algebraic expressions defined in the HMT model. If we were to use a more numerically complex mechanics model (see [15] for an example), we would expect an increase in computational time for the mechanics solution step similar to that shown in Figure 8 between the activation solution steps for the FK–HMT model and those for the N–LRd–HMT model. The time-dependent, dynamic component of active force development is included in the starting solution of each mechanics time step, but we neglect changes in the dynamic component during the application of numerical length perturbations in order to facilitate the numerical techniques being used to solve the finite elasticity model.

In summary, the relatively simple FK–HMT-based model of a 4 × 4 × 4-mm tissue cube required 550 MB of memory and 43 hours of compute time (wall clock) to perform a one-second contraction simulation. The N–LRd–HMT-based model, which includes a biophysically detailed model of cellular electrophysiology, required 1.5 GB of memory and 61 hours of compute time. These numbers obviously have a serious potential impact when simulations on more complex geometries within a reasonable time period are under consideration. We have, however, recently extended this work to a simplified model of the left ventricle (LV) of the heart [41].

Figures 9 and 10 present the results of a simulation of the contractile portion of the cardiac cycle in a simple left-ventricular geometry. At the cellular level, this model consists of the FK–HMT electromechanics model described above. The pole-zero constitutive law [47] is used to describe the passive mechanical behavior of the tissue microstructure. A constant-volume cavity constraint provides the dynamic pressure load applied to the endocardial surface of the LV model during the isovolumic contraction and ejection phases of the cycle. During ejection, the cavity model extends beyond the basal plane of the LV endocardial surface, corresponding to the quantity of blood exiting the ventricular cavity.

Figure 9 Figure 10

The model presented in Figures 9 and 10 consisted of only 320,000 cells, much fewer than would be required for a spatially converged solution of the electrical activation model, but adequate for this initial study. This simulation was performed on the same IBM p690 machine described above, and memory requirements peaked at 4 GB. The one-second simulation required 304 hours of compute time (wall clock) running on 12 processors.

From the comparison of the FK–HMT and N–LRd–HMT cube models above, we predict that models of the LV using the N–LRd–HMT model would require approximately 12 GB of memory and 430 hours of compute time. However, the more complex biophysical N–LRd–HMT model may be unsolvable on such a coarse discretization for the electrical activation model. As we move to more realistic geometrical models, the computational requirements would increase even further.

We have developed a novel computational modeling framework for the simulation of cardiac electromechanics. We have shown how our framework can be used for simulations that use geometric models that vary from tissue-block models to ventricular models. This framework enables the development and testing of new hypotheses associated with ventricular pacing, myocardial ischemia, and defibrillation. The most significant drawback in the framework highlighted by the simulations presented above is the sheer amount of computational time required to obtain these results. This computational limitation currently precludes the embedding of detailed cellular models in more anatomically based geometries. By combining the N–LRd–HMT model with existing anatomically based models of the heart [48, 49], a model could be created that accounts for mechano-electrical feedback and stretch dependence in an anatomically based ventricular geometry via the tight-coupling solution procedure outlined above. This coupling between cellular contraction and activation has been proposed as a possible mechanism underpinning the heterogeneous electromechanical delay that is required to simultaneously produce physiological spatio-temporal sequences of both activation and contraction [39].

Our framework is sufficiently general that, through the use of CellML, new cell models can be accommodated without changing the existing software tools. Thus, as the rapid increase of computational resources continues, we anticipate that we will soon be using this framework with more detailed anatomy and biophysically based cell models. As an illustration of this, we performed the above cube simulations using a new IBM pSeries* 595 with 256 GB of memory and 64 1.9-GHz IBM POWER5* processors. The N–LRd–HMT-based simulation, which required 61 hours of compute time on the p690 machine using eight processors, ran in eleven hours using 32 processors on the p595. This shows the increased computational performance of the newer machine and provides the platform to perform more complex computations. In addition to simply obtaining larger and faster computers, we are also investigating many areas of algorithmic and software design that will provide even greater improvements in computational cost. For example, speed increases may be obtained by compilation of CellML models into optimized descriptions using lookup tables and partial evaluation, the use of multigrid techniques and adaptive local mesh refinement, and altering code design to make use of distributed massively parallel processing environments.

Further application of the research described above has also begun in other organ systems that involve electromechanics. For example, the musculoskeletal system is a large system in which modeling of electromechanics is important to the understanding of function. While the detail of the skeletal muscle electromechanics differs from that of cardiac muscle, the underlying methods we have developed for the cardiac models are equally applicable to skeletal muscle [50]. With CellML, it is straightforward to replace the cardiac cellular models with skeletal muscle models, and simulations of knee flexion, for example, have been performed [51].

Like skeletal muscle, smooth muscle undergoes active contraction in response to electrical stimuli. While skeletal and cardiac muscle involve somewhat similar time scales, smooth muscle is significantly slower. Again, the underlying modeling framework developed in this work is capable of representing models of contracting smooth muscle. Initial investigation of this class of muscles is currently underway in the modeling of an active bronchial airway [52]. The computational modeling framework discussed in this paper provides a simulation environment for the integration of various physics approaches over multiple space and time scales. As such, the framework provides an ideal platform for the development of innovative models of human physiology in order to aid medical sciences in clinical diagnosis and drug discovery.

This work was supported by the Wellcome Trust, the Royal Society of New Zealand (Centre for Molecular Biodiscovery), and the Marsden Fund.

*Trademark, service mark, or registered trademark of International Business Machines Corporation.
**Trademark, service mark, or registered trademark of The MathWorks, Inc., Wolfram Research, Inc., or the Mozilla Foundation in the United States, other countries, or both.

Received September 27, 2005; accepted for publication December 20, 2005; Published online June 27, 2006.