## Abstract

We present the current state of the development of the SAPHIR project (a Systems Approach for PHysiological Integration of Renal, cardiac and respiratory function). The aim is to provide an open-source multi-resolution modelling environment that will permit, at a practical level, a plug-and-play construction of integrated systems models using lumped-parameter components at the organ/tissue level while also allowing focus on cellular- or molecular-level detailed sub-models embedded in the larger core model. Thus, an *in silico* exploration of gene-to-organ-to-organism scenarios will be possible, while keeping computation time manageable. As a first prototype implementation in this environment, we describe a core model of human physiology targeting the short- and long-term regulation of blood pressure, body fluids and homeostasis of the major solutes. In tandem with the development of the core models, the project involves database implementation and ontology development.

## 1. Introduction

The realization of the physiome project depends essentially on the development of a toolbox that will enable cooperation and collaboration among widely dispersed laboratories working on different platforms, in various modelling environments, on a variety of problems in different organs/tissues/cells at many different scales. The STEP1 Roadmap (see Fenner *et al*. 2008) identifies the items of such a toolbox, some pieces of which are already well advanced (as is evident from the articles in this special issue), others are in active development and still others are being planned or remain to be conceived. The proving ground of these tools will be the eventual constitution of a physiome modelling environment allowing the standardized reuse of existing models and their fusion into composite models via a plug-and-play (‘services’) approach.

One of the physiome toolbox proposals described in the STEP Roadmap is the development of ‘core models’. In the interest of providing a context for an *in silico* exploration of system-level integrated feedback and regulation loops, compatible with focused high-resolution detailed descriptions of local phenomena, these core models will consist of a *basic set* of low-resolution modular models of integrated function built around available experimental descriptions of global behaviour. In addition to the core set of low-resolution modules, provision will be made for the inclusion of higher resolution local sub-models that can be combined to customize the global model as replacements for one or more of the core modules. The core models will call on Virtual Physiological Human (VPH)/physiome toolbox utilities such as mark-up languages (CellML, FieldML: Cuellar *et al*. 2003; SBML: Hucka *et al*. 2003; or variants/extensions of these), thus serving not only as examples of their implementation but also as a testing ground. The goal is to keep the basic core models compact enough to ensure fast execution time and yet to allow elaborate detailed sub-models of target tissues or organs while maintaining the system-level regulatory compensations.

Here, we present the current state of the development of the Systems Approach for PHysiological Integration of Renal (SAPHIR) project2 (Thomas *et al*. 2007), a multi-resolution core modelling environment, with the implementation of a prototype core model of human physiology targeting the short- and long-term regulation of blood pressure, body fluids and homeostasis of the major solutes. Also included are the main regulatory sensors (baro- and chemoreceptors) and nervous (autonomic control) and hormonal regulators (antidiuretic hormone, aldosterone and angiotensin). The resulting modelling resource will be provided as an open-source project available to the general community as part of the physiome effort.

While models have previously been developed with a similar systems approach (Srinivasan *et al*. 1996; Abram *et al*. 2007), our model aims to go further towards the provision of an open-source, modular resource, amenable to the plug-and-play exploration of the systems implications of detailed and focused models of individual subsystems. We intend to produce scaled versions of the basic core model for human, of course, and also for rat, mouse and dog, since these are the principal experimental models.

Rather than starting from scratch, two legacy integrated models have served as the starting point, namely the classic model of Guyton *et al*. (1972), which focused on blood pressure regulation, and the model of Ikeda *et al*. (1979), based on Guyton's models but extended to focus on the overall regulation of body fluids. We have functioning re-implementations of the Guyton *et al*. and Ikeda *et al*. models in Fortran, C++, Simulink and also using the ODE solver Berkeley Madonna. The set of parameter values, along with related work from the literature, is being collected in a quantitative systems database (QSDB) that is an extension of our generic QxDB/quantitative kidney database (QKDB) database environment (Dzodic *et al*. 2004; Ribba *et al*. 2006). This also involves ontology development.

## 2. Model descriptions

### (a) The 1972 Guyton model of overall circulation (G72)

Although it was published over 30 years ago, the original model of Guyton *et al*. (1972) remains a landmark achievement. Furthermore, with the rise in the last 10 years of systems physiology, the Guyton model has attracted renewed attention (see Simanonok *et al*. 1994; Wabel & Leonhardt 1998; Nguyen *et al*. 2008). This model was the first ‘whole-body’, integrated mathematical model of a physiological system. It allows for the dynamic simulation of systemic circulation, arterial pressure and body fluid regulation, including short- and long-term regulations. From a modelling standpoint, it is actually a composite approach, since the model uses both ‘exact’ physical and physiological laws (explicit and validated) and ‘curve fits’ of experimental data or simply tabulated data, e.g. left ventricular output as a function of systemic arterial pressure, LVM=*f*(PA), is given as a piecewise linear graph.

The 1972 Guyton model consists of 18 modules (350+ elementary blocks), containing approximately 160 variables, including more than 40 state variables. The model contains a total of approximately 500 numerical entities (model variables, parameters and constants). A current version used for teaching and equipped with a sophisticated graphical user interface, developed continuously by members of the original Guyton laboratory, now includes approximately 4000 variables, including many physiological entities or processes that were as yet unidentified in 1972 (Abram *et al*. 2007). However, this more elaborate model is less well suited to our purposes than the 1972 model, since there is as yet no complete description of its formulation.

In essence, Guyton's original model is constructed around a ‘central’ *circulatory dynamics* module in interaction with 17 ‘peripheral’ modules corresponding to physiological functions (such as *pulmonary dynamics and fluids* or *non-muscle oxygen delivery*; see Guyton *et al*. 1972). An examination of the original code or published diagram reveals that, in addition to its interconnected module structure, the model is characterized by a wide range of time scales in the different modules, ranging from 5×10^{−4} min (*autonomic control*) to 10^{4} min (*heart hypertrophy or deterioration*). Moreover, the original Fortran program code contained several programming and computation-related features that required special treatment.

The basic conservation relations of the 1972 Guyton model concern sodium, potassium, oxygen, protein, water and the blood volumes in the circulatory loop. These are given by 20 ordinary differential equations, which we summarize here to give an idea of how the model was developed, without going into detail for each of the balance relations. Total extracellular fluid sodium (NAE) depends on the balance between dietary sodium intake (NID) and urinary sodium excretion (NOD; where the prime (′), here and below, indicates the time derivative),(2.1)Extracellular (KE) and intracellular (KI) potassium levels depend on the rates of potassium intake (KID) and renal loss (KOD) and on the rate of transfer from extracellular to intracellular compartments (KCD),(2.2)(2.3)Muscle cell O_{2} (QOM) depends on the rate of O_{2} delivery (RMO) and its rate of usage (MMO). Similarly, for non-muscle cell O_{2} (QO2), where DOB and MO2 are rates of delivery and usage in these cells,(2.4)

(2.5)

As shown in figure 1, protein is distributed among plasma (PRP), interstitial fluid (IFP), interstitial gel (GPR) and pulmonary fluid (PPR). The variables PPN, PPO, DPC, DPL and GPD represent the rates of exchange between compartments, as shown in figure 1, and DLP and DPO are the rates of protein formation by the liver and total rate of loss of protein, respectively,(2.6)(2.7)(2.8)(2.9)Figure 2 shows the distribution of body water among: blood plasma (VP) and red cells (VRC), which together make up total blood volume (VB), pulmonary free fluid (VPF), interstitial gel (VG), interstitial free fluid (VIF) and cell volume (VIC). The variables PFI, PLF, VTC, VTL, VGD, VID, RC1, and RC2 represent the rates of exchange between compartments, as shown in figure 2, and TVD and VUD are the rates of drinking and urinary volume excretion, respectively,(2.10)(2.11)(2.12)(2.13)(2.14)(2.15)Figure 3 shows the distribution of blood as it flows through the main compartments of the general circulation, namely right atrium (VRA), left atrium (VLA), systemic arteries (VAS) and veins (VVS), and pulmonary arteries (VPA). The variables QVO, QRO, QPO, and QLO represent blood flow at various points along the circulation. BFM and BFN are the muscle and non-muscle blood flow, respectively, and RBF is the renal blood flow. The terms PLA, PPA, PRA, PA and PVS represent the five compartmental pressures, relative to atmospheric pressure. In the model, the rates of volume change of these main compartments are given by(2.16)(2.17)(2.18)(2.19)(2.20)where the last term in each equation represents the redistribution of changes of total blood volume (VB; due, for example, to haematopoiesis or, on the contrary, to haemorrhage) among the various compartments (the fractions sum to 1.0).

### (b) Ikeda *et al*. (1979) model of body fluid and acid–base regulation

Ikeda *et al*. (1979) published ‘A model of overall regulation of body fluids’, which was an extension of the 1972 Guyton model focused on acid–base regulation. Their model includes subsystems for the circulation, respiration, body fluids, renal function, and acid–base balance. This model is well adapted for problems concerning body fluid disturbances and fluid therapy. The model consists of seven blocks (figure 4) grouped into three categories, which are as follows.

The ‘circulation and body fluids’ category contains blocks 1, 3 and 4, which, respectively, correspond to the cardiovascular system, the extracellular space, and the intracellular space and electrolytes.

In the cardiovascular system (

*block*1), the steady-state values of cardiac output and of the mean pressures in systemic and pulmonary arteries and veins are calculated as linear functions of the circulatory blood volume, on the basis of the Starling mechanism's linear approximation. The model is concerned only with normal states of the heart, and the circulatory subsystem is simplified when compared with Guyton's 1972 model—rapid transients such as autoregulation or stress relaxation are ignored.*Block*3 deals with water intake and output and with the distribution of extracellular fluid. Protein flows in plasma and interstitial fluid are also represented.*Block*4, related to intracellular space, deals with the dynamics of the intracellular fluid volume and with the intracellular acid–base balance. The concentrations of sodium, potassium, glucose, urea, and mannitol are calculated. This block is more detailed than in Guyton's model; indeed, Ikeda*et al*. wanted a model that would reproduce/simulate the findings in experimental metabolic acidosis and alkalosis.The ‘respiration’ category consists only of the respiratory system (

*block*2). In this block, the alveolar ventilation is calculated as a function of pH and partial pressures of CO_{2}and O_{2}in the blood. These partial pressures in the arterial blood as well as in the alveoli and the concentrations of CO_{2}and O_{2}in the arterial and venous blood are evaluated in this block thanks to dynamic relations.Renal function, which was relatively simple in Guyton's 1972 model, is more detailed here. This category is composed of three blocks:

*blocks*5 and 6, related to the kidney, and*block*7, described as the controller of the renal function.*Block*5 deals with the renal excretion rate for bicarbonate, calcium, magnesium, phosphate, sulphate, chloride, and organic acids. Reabsorption and excretion of water, sodium, potassium, ammonia, and titratable acid are calculated in*block*6. Urine pH is evaluated from the excretion rates of phosphate and organic acid according to the Henderson–Hasselbalch equation.*Block*7 gives the glomerular filtration rate and the control signals that operate on the kidney, such as antidiuretic hormone or aldosterone. The renin–angiotensin–aldosterone system is not detailed.

For each block, inputs either come from the other blocks of the model or are fixed or constant parameters. Concerning the outputs, some affect the other blocks of the model and some are independent of the other blocks. According to the authors, this partition was made to minimize the interaction between the blocks.

### (c) Development of detailed sub-models

In parallel with the development of the core model, we are building a collection of detailed sub-models for the major organ systems involved in blood pressure regulation and systems-level fluid and solute homeostasis. This involves the adaptation of legacy, published models at many scales and at various levels of resolution. The adaptation consists of implementation in the ‘language’ of the SAPHIR toolbox, with a standardized specification of input and output variables, including the appropriate handling of boundary conditions in the case of coupled multi-resolution sub-models. Several such extensions are under development for the various organs, but here we mention only the renal sub-models.

#### (i) Renal sub-models

The kidney is the main actor in the long-term regulation of blood pressure and hydroelectrolytic equilibrium (Guyton 1990*a*; Cowley 1997; Ferrandi *et al*. 1999; Bianchi *et al*. 2005). As such, it plays a key role in all forms of hypertension (e.g. Lindner *et al*. 1992; Alvarez-Guerra & Garay 2002; Hannaert *et al*. 2002) and, in particular, in sodium-dependent forms of hypertension (Campese & Park 2006). On the other hand, hypertension accelerates diseased kidney function loss, and antihypertensive treatment slows renal disease progression.

The primary role of the kidney module (KM) in the core model is to represent the kidney-dependent functionality related to the regulation of blood pressure, volaemia and sodium balance, i.e. to represent the renal function curve, also called the pressure–natriuresis curve (Guyton 1990*b*; Hall *et al*. 1996). Detailed sub-models for the myriad of local mechanisms at the level of glomerular filtration, tubuloglomerular feedback, transport along individual nephron segments, cortical and medullary microcirculation, the role of the medulla in urine concentration, etc. can be implemented and used to replace part or all of the basic KM in order to address specific problems. These detailed local models may of course be based on legacy models (Thomas *et al*. 2006*a*,*b*) or developed from scratch.

## 3. Modelling and knowledge management infrastructure

Here, we describe the modelling and simulation tools and our parameter database, and we discuss the ontology tools we are developing as a support for these.

As stated previously, a core model consists of a number of related modules, each having clearly specified input and output variables. Within each module, the output variables are calculated based on the values of the input variables, either by some simple transfer function or by a more-or-less detailed model of the subsystem involved. The sub-models may be developed in any of a number of different formalisms, e.g. differential algebraic equations, multi-agent systems, cellular automata, etc. Consistent with the object-oriented spirit of the core modelling environment, the only restriction on each sub-model is that it must respect the imposed input/output characteristics in order to correctly interface with the rest of the systems model. Within a given sub-model, ‘local’ variables may be defined, permitting to follow, during actual *in silico* simulations, any additional appropriate processes taking place within the considered subsystem. For example, a hybrid (multi-formalism) electrophysiological model of an ischaemic cardiac tissue can be constructed as a set of coupled sub-models representing healthy and ischaemic cells (Defontaine *et al*. 2004, 2005). Healthy cells can be represented by means of a simplified (discrete) automaton model, while ischaemic cells can be represented with a detailed (continuous) model. We can choose to store as output only the state variables of some ischaemic cells (local variables of the continuous sub-models) and the global activity at the tissue level (the EGM potential generated by the set of cells).

This approach requires a robust and flexible simulator architecture in order to keep the many processes in step with each other while minimizing computation time. Numerical techniques must also be appropriate for the *stiff* nature of the problem, since some processes are orders of magnitude faster than others in this system (e.g. in the above example, a fourth-order Runge–Kutta (RK) integration method with an adaptive step size is used for detailed models, while a discrete-time simulator is used for simplified models). The LTSI team (Rennes) has developed such an environment, called Multi-formalism Modelling and Simulation Library (M2SL), based on the co-simulation principle (Vangheluwe *et al*. 2002) and has demonstrated its use for several problems in cardiovascular physiopathology and epileptology (Hernández *et al*. 2002; Le Rolle *et al*. 2005; Wendling *et al*. 2005). These techniques scale up nicely to the systems model level of the present project.

### (a) M2SL

M2SL is a collection of generic C++ classes that has been specifically developed for modelling and simulating complex systems using a combination of different modelling formalisms. The library is based on the work by Zeigler *et al*. (2000) and Vangheluwe *et al*. (2002), who proposed a distributed co-simulation approach to couple multi-formalism models. In this approach, each model object is composed of a hierarchy of coupled sub-models that may be defined with different formalisms (discrete, continuous, etc.) or different resolutions (sub-cellular, cellular, tissue, organ, or system levels). A parallel object hierarchy, in which a specific local simulator is coupled to a corresponding model object, is defined automatically at run-time: a continuous (sub)model will be associated with a continuous local simulator, and a discrete (sub)model will be associated with a discrete local simulator. Local simulators will adapt their properties (such as the integration step size in the continuous case) according to the dynamics of the corresponding model.

The coupling between local simulators is performed on global simulators or ‘coordinators’. The main difficulty associated with the development of such coordinator objects concerns the definition of completely general criteria for spatial coupling and temporal synchronization. Depending on the formalisms involved, spatial coupling can be done by means of sampling-and-hold, interpolation (linear, splines, etc.), or quantization methods. Three types of temporal synchronization have been implemented in the library: (i) simulation and global synchronization with a fixed step size for all elements, (ii) global synchronization at a fixed time step and local simulation with an adaptive step size, and (iii) local simulation with the smallest time step and global synchronization with a variable step size. The co-simulation approach presents several advantages with respect to classic modelling and simulation systems: (i) it is particularly well suited for the creation of generic multi-formalism modelling tools, (ii) it facilitates distributed execution in a computer cluster, and (iii) it eases the extension of the model to a multi-resolution representation.

### (b) Parameter database

Given the need for reliable databases of experimentally determined parameters in order to underpin and streamline the development of new models, we have developed a QSDB, an extension of our QKDB (Dzodic *et al*. 2004; http://physiome.ibisc.fr/qkdb/), to include the other relevant organ systems. These are Apache/PHP/MySQL databases. The content is currently more than 8000 entries from approximately 300 articles, but this is still far from complete coverage of the relevant data. The ultimate success of this effort as a central source for quantitative physiological data will depend on participation (i.e. data entry) by a wide cross section of the physiology community. Since the data model of QKDB was designed specifically to allow flexible extension of data descriptors, we also extracted a generic version, QxDB (Ribba *et al*. 2006), that can serve as the basis for similar databases in new domains.

### (c) Ontologies

An ontology, in computer science, is a hierarchical representation of a set of concepts within a domain and the relationships linking those concepts. It provides the semantic underpinning necessary for the automated treatment of data or, in our case, mathematical models and their parameters. A number of medical and biological ontologies are available through the OBO foundry website (http://www.obofoundry.org) that has links to, for instance, the Gene Ontology (a set of three ontologies that aim to gather all the terms related to molecular functions, biological processes, and cellular components) and the National Library of Medicine's Unified Medical Language System, actually a metathesaurus. Specifically, biomedical ontologies are the subject of a recently created National Centre for Biomedical Ontology (Rubin *et al*. 2006). A helpful discussion of ontologies in the physiome can be found in Crampin *et al*. (2004).

Within our project, ontologies serve several needs. First, the entries in QKDB require unambiguous descriptors. For some categories of terms, such as anatomy, there exist ontologies that can be used, such as the Foundational Model of Anatomy for anatomical terms (with extensions for missing details). For certain other domains, however, there is, to our knowledge, no existing ontology—this is the case especially for terms describing physiological parameters, such as permeability, flow rates, viscosity, channel density, etc. Furthermore, among workers in different communities, a given term may refer to different concepts or, more frequently, there are several synonyms for the same concept or object. In our project, we need to use different ontologies that currently have no direct links but that serve related domains. In order to use the needed ontologies in an integrated way, ‘bridges’, i.e. mappings, between ontologies must be built by defining the semantic relations between concepts in the different ontologies. Mappings are often established manually by experts, but, owing to the size and the number of ontologies, there is a need for automatic (or semiautomatic) mapping. Several algorithms and tools have been proposed in the literature (McGuinness *et al*. 2000; Noy & Musen 2004; Ehrig & Sure 2005; Lambrix & Tan 2006; Bouzeghoub *et al*. 2008), but due to the complexity of the task, the problem is still open. We are developing a system based on a multi-agent cooperative approach and on machine learning methods in order to combine efficiently several similarity measure methods that exploit the syntactic, semantic, and structural information contained in ontologies for related domains. This method takes into account each characteristic of the concepts and the relationships defined by the different ontologies. External static dictionaries, dynamic dictionaries, and distance measures (metrics based on ontology tree branchings) are used to infer the connections that can be built among these ontologies. Of course, false positives will require curation by biological experts, but this will be less time and energy intensive than a complete manual resolution of the problem.

Second, concentrating on the domain of blood pressure regulation and fluids homeostasis, we have built an ontology based on all the variables and parameters in the Guyton & Ikeda models. Although certainly not exhaustive for the domain as a whole, this approach at least includes a self-consistent set of terms covering the main concepts, as represented by these validated legacy models. As detailed, local sub-models are added to the core model collection, the additional concepts will be added. The terms of this ontology are used to tag the variable names in the core collection of sub-models, as well as the equations in which they are used. This will also serve as a link to entries in our database relevant to the model variables and parameters. The ontology thus serves as a key tool for the navigation of the core models, as well as being an extensible description of the field and a valuable complement to QKDB/QSDB.

## 4. Core model results

We have implemented the Guyton model in two modular environments, M2SL and *Simulink*; the Ikeda model is implemented under *Berkeley Madonna*. All three environments allow straightforward replacement or revision of model details. Preference will depend on users' backgrounds.

### (a) G72 under M2SL and Simulink

#### (i) G72 under M2SL

In order to implement the 1972 Guyton model (G72) using M2SL, one model object was defined for each of the main ‘modules’. A class diagram of the implementation is presented in figure 5. The class ‘Guyton’ is the main model class, and it creates at least one instance of all other classes as sub-models. As explained previously, a continuous simulator object is created for and associated with each of the model objects, and the Guyton object is coupled to a coordinator simulator that will handle input–output and temporal coupling between all sub-models. The adaptation properties of each simulator facilitate the implementation of the model and minimize simulation time. This has been particularly useful for the implementation of the G72 model, since various blocks operate with different time constants. For example, blocks representing haemodynamic function are characterized by shorter time constants than those representing hormonal regulation. In order to take into account this heterogeneity, the original Fortran implementation manually imposed a shorter integration step size for the haemodynamic and the autonomic regulation blocks than for the other parts of the model. In the M2SL implementation, the simulators associated with the haemodynamic and autonomic nervous system regulation blocks automatically use an integration step size approximately 100 times shorter than the other block simulators.

In order to validate the M2SL C++ version of the model, we simulated two *in silico* experiments described in the 1972 Guyton *et al*. paper and compared the results with the output from the original Fortran program in Guyton's laboratory (provided, along with the original Fortran code, by R.J.W., who worked in Guyton's laboratory at the time). The first, which we call benchmark 1 (BM1), is the simulation of hypertension in a salt-loaded, renal-deficient patient by first reducing the renal mass by 70% and later increasing the salt intake fivefold at, respectively, 2 hours and 4 days after the beginning of the simulation. The simulations correctly predict that the reduction in renal mass induces a decrease in cardiac output and an increase in peripheral resistance and arterial pressure. In response to the increased salt load, the extracellular volume, blood volume and cardiac output rise while the total peripheral resistance falls (figure 6).

The second *in silico* experiment, benchmark 2 (BM2), is the simulation of sudden severe muscle exercise and takes place over a much shorter time scale than BM1 (9 min instead of 8 days). The exercise is simulated as follows: at *t*=30 s, the exercise parameter was changed to 60 times its normal value (EXC=60.0), corresponding to an approximately 15-fold increase in the whole-body metabolic rate (for this step, the time constant for the local vascular response to metabolic activity was reduced by 1/40 (A4K=0.025) and certain numerical damping factors were modified); then, at *t*=2 min, the value of *EXC* was reset to normal (EXC=1.0). As can be seen in figure 7, within seconds of the onset of exercise, cardiac output and muscle blood flow increased considerably. The urinary output fell to its minimum level while arterial pressure rose moderately. The muscle cell and venous PO_{2} fell rapidly. The muscle metabolic activity showed an instantaneous increase, but then decreased considerably owing to the development of a metabolic deficit in the muscles. When the exercise was stopped, the muscle metabolic activity fell below normal, but cardiac output, muscle blood flow and arterial pressure remained elevated for a while as the person was repaying his/her oxygen debt.

Figures 6 and 7 show the close match between M2SL simulations and the results from the original Guyton model.

#### (ii) G72 under Simulink

*Simulink*, a block-based language for describing dynamical systems, was also used as a modelling and simulation platform (version 6.1, integrated with Matlab, Release 14, The MathWorks, Nattick, MA, USA). It is an interactive and graphical environment dedicated to the multi-domain simulation of hybrid continuous/discrete systems. During simulations, model and block parameters can be modified, and signals can be easily accessed and monitored. In the models described below, numerical integration was performed using ‘ode15S’ (a Matlab library) with a variable step size (maximum step size, *auto*; relative tolerance, 10^{−3}).

As for the M2SL implementation, we based our Simulink version of G72 on the detailed control system diagram published in the original article and on the Fortran code of the original model.

First3, code operations and routines from the computer program were rendered into the Simulink graphical description, i.e. elementary blocks and subsystems connected by appropriate signals. This led to a structure similar but not identical to that of the original model diagram, i.e. one central *circulatory dynamics* module that shares information (via common global variables) with peripheral subsystems, such as capillary membrane dynamics, kidney function, or hormonal control.

Second, subsystems were treated as ‘atomic subunits’ (by selecting the proper option in the ‘Subsystem parameter’ dialogue box). This causes Simulink's solver to treat each subsystem as a computation unit, i.e. to perform all calculations within that subsystem before going to the next (otherwise, the execution of blocks in a given subsystem could be interleaved with the execution of external blocks). Technically, this led us to use the ‘Discrete integrator’ block. This Simulink predefined block allows explicit specification of its ‘sample time’, which enabled us to use the state variable integration scheme of the original Fortran program. However, this functional structuration still did not result in a functioning Simulink model, i.e. we observed a lack of convergence due to oscillations and other run-time errors.

Third, the execution order of the original Fortran program was imposed, using the option proposed by Simulink to prioritize the execution order. This was done on both subsystems and their lower level algebraic groups of blocks. This modification greatly stabilized the model and allowed it to run to completion, under the ‘baseline’ settings. A quasi-steady state was reached within minutes of simulated time (stable for at least 2×10^{4} min simulated time; time derivative of state variables in the range 10^{−12} to 10^{−6}).

Fourth, in order to deal with the marked stiffness of the model, subsystems were regrouped into ‘fast’ blocks, with characteristic time constants in the 10^{−4} to 1 min range, and ‘slow’ blocks (1–10^{4} min). The modules *circulatory dynamics*, which computes pressure, volume, and flow in the various circulatory compartments, and *autonomic control*, which computes regulatory influences by the autonomic system, are in the Fast subsystems category. All other modules from Guyton's model are in the Slow category; principally, these 21 sub-models treat extravascular compartments, fluid movements and regulations, pulmonary dynamics, oxygen delivery, and blood flow control, kidney function and hormonal control (antidiuretic hormone, angiotensin and aldosterone).

Numerical integration was performed in Simulink using (i) ‘discrete’, ‘variable step’ and ‘max step size=0.1 min’ as solver options and (ii) ‘forward Euler’ method for discrete integrator blocks within the model.

As with the M2SL model described above, we used the output from Guyton's *in silico* experiments (Guyton *et al*. 1972) as benchmarks to validate the implementation of Guyton's 1972 circulatory model (‘G72’) with Simulink.

Benchmark results were similar to those obtained with M2SL and are not shown. For BM1, a quantitative comparison of Simulink versus benchmark simulations yielded an average value of relative r.m.s.e. (*r*r.m.s.e.) of 0.050±0.079 (mean SD, *n*=9 variables). For BM2 (severe muscle exercise), *r*r.m.s.e. was 0.049±0.061 (mean SD, *n*=8 variables).

### (b) Ikeda model under Berkeley Madonna

The model of Ikeda *et al*. (1979) was very clearly and completely described in the original paper. The model was illustrated with one block diagram (figure 4), defining the global structure and the links between the inputs and outputs of the blocks. Each block of the model is also specified by one schema (an example corresponding to *block* 3 is shown in figure 4*b*) composed of symbols that are mainly mathematical operators (addition, product, and integration). Besides those symbols used for drawing the model, several black boxes represent certain main functions, whose mathematical expressions are also given in the article's appendix.

The entire Ikeda model is defined by more than 200 variables and constants and consists of approximately 130 nonlinear differential or algebraic equations. Some of these equations are given explicitly in the main body of the article. Two examples from *block* 3 are as follows. Equation (4.1) calculates the pressure of the interstitial fluid (PIF) starting from its volume (VIF),(4.1)where *x*=VIF/VIF_{0}. Equation (4.2) calculates the rate of the lymph flow (QLF) according to the VIF,(4.2)

Other equations are obtained thanks to the various mathematical symbols in the block schemas. Two examples extracted from *block* 3 are as follows. Equation (4.3) gives the VIF from the capillary filtration rate (QCFR), the rate of the lymph flow (QLF) and the rate of the water flow into the intracellular space (QIC),(4.3)Equation (4.4) evaluates the volume of the extracellular fluid (VEC) as the sum of the plasma volume (VP) and the interstitial fluid volume (VIF),(4.4)

The parameter values and normal values of variables detailed in the original article were determined by Ikeda *et al*. using experimental results from the literature and supplementary experiments. The model assumes a healthy human male of approximately 55 kg in weight.

The Ikeda model was initially written in the Fortran IV language. We have implemented a new version using the *Berkeley Madonna* software (Macey & Oster; http://www.berkeleymadonna.com), a fast and convenient differential equation solver with a simple programming language. Of the several integration methods proposed, the fourth-order RK method worked well for the computation of the Ikeda model. Typically, the computation time was a few seconds for a few hours of simulated experiment.

The complete transposition of the model was successful after the correction of several bugs and mistakes detected in the original paper (continuity problems in mathematical expressions, inconsistencies between formulae explicitly written in the article and symbols in the blocks schemes or between the normal values displayed on the blocks and those listed in the appendix).

A few lines of the Berkeley Madonna implementation for *block* 3 of the Ikeda model are presented below, corresponding to the four differential and algebraic equations described above.

#### (i) Simulations with the Ikeda model

Once implemented under Berkeley Madonna, we repeated the simulations described in the original paper. In Ikeda *et al*. (1979), these were set against the data from the literature to show that the model behaviour agreed well with the experimental results pertaining to body fluid and acid–base disorders.

Here, we present one simulation: the respiratory responses to increased CO_{2} concentration in inspired air. We simulate, during 1 hour, the transient response of the respiratory parameters to the inhalation of 5% CO_{2} in air over 30 min (FCOI=0.05 from *t*=5 to 35 min). We observe the alveolar ventilation (VI), the pressure of CO_{2} and O_{2} in the alveoli (PCOA and PO_{2}A), and the concentration of bicarbonate of the extracellular fluid (XCO_{3}). The results are compared with those in the figures of the original article (figure 8). The black lines correspond to our Berkeley Madonna simulations and the grey dashed lines have been extracted from the original paper. The results are quite similar in both simulations; the *r*r.m.s.e. computed from the simulated and original signals and averaged over the four variables is equal to 0.0305. Similarly, a good agreement was obtained with simulations of the other *in silico* experiments presented in Ikeda *et al*. (1979; results not shown).

## 5. Discussion and conclusions

In keeping with the spirit of the physiome, the SAPHIR core modelling environment is a step towards a collaborative platform for the construction of open-source, multi-scale, multi-mode models of integrated physiological function in a reasonably low-CPU computing environment compatible with future applications in a clinical setting. A framework with similar goals was developed (in the UCL Beacon Project) and applied to hepatocyte glycogenolysis (Hetherington *et al*. 2007).

This approach to the VPH is, in a certain sense, the antithesis of the chimeric idea of an all-inclusive model that would supposedly include everything: gene expression and protein folding; metabolism; and detailed models of all 200+ cell types, of all tissues and organs, of the immune system, and, even by extension, of the whole nervous system and the brain. As argued forcefully by Rosen (2000), such an all-encompassing, ‘largest model’ of a *complex system* such as the human body is impossible even in principle. Rather, we must recognize that only an infinity of local models can approach a complete description of a complex system, much as a collection of local planar maps approaches a description of a spherical world, with which one must shift from one local map to another to see details, while the global sphere maintains the overall context (Thomas 2007).

The SAPHIR core modelling environment is thus a multi-scale approach that will, when fully implemented, allow arbitrarily detailed local models at any scale, when implanted in a larger core model, to take advantage of full accounting for system-level feedback regulations when called for by the problem of interest. Certainly, one such problem is the regulation of arterial blood pressure, as was amply established by Guyton and his many collaborators, since their quantitative systems models led them to profoundly reorient our understanding of the causes of hypertension. This was our rationale for adopting the early Guyton *et al*. (1972) model and its extension to acid–base regulation by Ikeda *et al*. (1979) as the initial demonstrators for the core modelling environment. Obviously, ours is neither the first nor the last (re-)implementation of the Guyton models, which were not only implemented with a generic user interface (called MODSIM; see Montani *et al*. 1989), but have also been usefully deployed in a number of contexts, including research on the physiological consequences of weightlessness in manned space flight (White *et al*. 1991, 2003). An implementation of the Guyton model from the published diagram (Guyton *et al*. 1972) with Simulink has recently been published (Kofranek & Rusz 2007), not to mention the current very elaborate and comprehensive Quantitative Human Physiology model at the University of Mississippi, directly descended from the early Guyton models by the colleagues of Guyton (Abram *et al*. 2007).

The originality of our core model implementation is our commitment to provide documentation for each basic module, interactive modification of any aspect of the model parameters or equations and, especially, the possibility of contributing new detailed local sub-modules to complement or replace selected ‘basic modules’.

Crucial to this approach is an underlying computational toolbox that is sufficiently robust to handle the wide range of spatial and temporal scales, and flexible enough to accept sub-modules in a variety of formalisms, including not only ODE/PDE dynamical systems models but also discrete modelling methods such as cellular automata. The M2SL toolbox was designed specifically to answer these criteria in the context of cardiac modelling and is currently our best choice.

Essential for the perennity of the environment are ontologies for the resolution of nomenclature conflicts, consensus on adequate mark-up languages for model description and data exchange, and curated databases to centralize the voluminous and heterogeneous data that underlie all such modelling efforts. Crucial also will be the establishment of procedures for model verification (mathematical) and validation (against experimental and clinical data). At present, these goals are only partially met, but they are the goals of the general VPH/physiome project on a European and international level.

## Acknowledgments

The authors thank Ricardo Garay for his help and expertise in hypertension. We thank Nadia Abchiche, Amel Bouzeghoub, Abdeltif Elbyed, Tarek Melliti, and Farida Zehraoui for their collaboration with the ongoing development of the ontology tools and algorithms.

This work was supported by the CNRS and the INSERM. S.R.T. gratefully acknowledges financial support from the French National Research Agency (ANR grants ANR-06-BYOS-0007-01 (SAPHIR) and ANR-05-BLAN-0247-06 (CorBioMath)), the European Union (IST-2004-IST-4 STEP), the French Ministry of Economy, Finance and Industry (Pôle System@tic, FAME2 project) and the Council of Essonne Region (Pôle System@tic, POPS project). R.J.W. would like to acknowledge partial support from Cooperative Agreement NNJ06HG25A between the Universities Space Research Association Division of Space Life Sciences and the National Aeronautics and Space Administration. This work was also aided by discussions within the GdR STIC-Santé CNRS 2647—INSERM.

## Footnotes

One contribution of 12 to a Theme Issue ‘The virtual physiological human: building a framework for computational biomedicine I’.

↵STEP: STrategy for the EuroPhysiome FP6-2004-IST-4.

↵ANR-06-BYOS-0007-01, http://saphir.ibisc.fr.

↵N.B. In the following, we use Simulink terminology: ‘block’ for elementary (pre-defined) graphical blocks (e.g. integrator) and ‘subsystem’ for groups of blocks (and/or subsystems) performing higher level calculations (e.g. logical test).

- © 2008 The Royal Society