1 Introduction

Computational modelling in physiology has contributed to many significant breakthroughs, but the models themselves have usually not become working tools for experimentalists nor even for other modellers outside the developer’s own group. We provide here a practical implementation of one of the classic and most complete models of body fluid and acid–base regulation, and we give several examples of the use of the model. We give the complete model description in the language of Berkeley Madonna, which is very easy to read and can readily be converted for other numerical solvers. Physiologists and clinicians will find this model easy to use, and this complete example will facilitate extensions in order to simulate related clinical situations or new experimental findings.

Inspired by the classic model of blood pressure regulation by Guyton et al. (1972a), Ikeda et al. (1979) adopted the same symbolic representation to illustrate model structure, but since their focus was on body fluids and acid–base balance, which have a slower time course than, say, autonomic regulation of cardiovascular variables, they simplified the representation of the cardiovascular system but greatly extended the renal and respiratory systems. Their model consists of a set of nonlinear differential and algebraic equations with more than 200 variables and has subsystems for circulation, respiration, renal function, and intra- and extra-cellular fluid spaces.

2 Materials and Methods

2.1 Model Description

The original article of Ikeda et al. (1979) describes the details of the model, so we will not give a complete description here (the program code, Online Resource 01, given in the Electronic Supplementary Material and described in the Appendix, has all the explicit equations); our implementation closely follows the description in their article, especially in their diagrams of the seven blocks that constitute the model, namely, the circulation and body fluids (blocks 1, 3, and 4), respiration (block 2), and renal function (blocks 5, 6, and 7). Initial values and many other details are given not only in the text but also on the diagrams and in the tables of the original article. Here, we give just a brief explanation of the basic content of the model and Ikeda et al.’s general strategy.

As in Ikeda et al. (1979), the model assumes a healthy male of approximately 55 kg body weight, and parameter values used here are those given in the original article. Calibration of the model for other body weights or for females would be a valuable extension of the model but is beyond the goals of the present work. Such extension would involve re-calibration not only of extracellular and intracellular fluid volumes (and thus with impact on solute contents of those compartments), but also of less straightforward parameters such as metabolic rate, respiratory volume, cardiac output, and the like.

The cardiovascular/circulatory (CV) system, quite complex in Guyton’s model, was considerably simplified by Ikeda et al. (1979) to a simple steady state that represents the system’s state after settling from transient local autoregulation or stress relaxation.

By contrast with the simplified CV system, and in keeping with their focus on acid–base and fluid physiology, Ikeda et al. (1979) included much more elaborate representations of the respiratory system, intracellular and extracellular electrolytes and solutes, and of course the kidney. For example:

  • Alveolar ventilation (VI) is calculated as a function of blood pH, \(\hbox {P}_{\hbox{CO}_{2}}\), and \(\hbox {P}_{\hbox{O}_{2}}\);

  • The blood buffer system is treated using the Henderson–Hasselbalch equation, an equation for the oxygen saturation curve, and an equation for the in vivo \(\hbox{CO}_{2}\) dissociation curve, thus the model takes account of the haemoglobin buffer system, the Bohr effect, and the Haldane effect;

  • The model treats intra- and extra-cellular electrolytes and acid–base balance and also glucose metabolism and insulin secretion—disorders of glucose metabolism can be modelled by varying the parameters CGL1, CGL2 and CGL3;

  • Plasma osmolality in the model depends on the concentrations not only of sodium, potassium, glucose, and urea, but also of mannitol, included in the model because of its frequent therapeutic use;

  • The renal blocks treat reabsorption and excretion not only of water, sodium, and potassium, but also of bicarbonate, calcium, magnesium, phosphate, and organic acids; proximal tubule reabsorption depends on volume expansion or pressure diuresis (THDF); aldosterone is assumed to act on the distal tubule to increase sodium reabsorption, decrease potassium secretion, and increase excretion of titratable acid; urine pH and excretion of ammonia, titratable acid, phosphate, and organic acids are included in the model; glomerular filtration rate (GFR), represented as a sigmoid function of arterial pressure, is controlled by extracellular volume (VEC) and depends on antidiuretic hormone (ADH) and aldosterone (ALD) and on THDF;

  • The renin–angiotensin–aldosterone system (RAAS) is represented here simply as a transfer function by which ALD secretion depends on extracellular fluid (ECF) potassium concentration, tubular sodium concentration, arterial pressure, and volume receptor signals.

In addition to this incomplete list, the model contains many other interesting features that the reader should glean from the original Ikeda et al. (1979) article.

2.1.1 Berkeley Madonna Description

Berkeley Madonna is a fast, robust, multi-platform solver of systems of ordinary differential-algebraic equations. Compared with other such solvers, it is extremely easy to program (a simple list of the equations in any order), has a very effective user interface for plotting or tabulating the results and varying the parameters (simple “sliders” can be easily defined to vary individual model variables or parameters, with instant re-run of the model), and it has proven to be very fast compared to other solvers we have used.

3 Results

To demonstrate several interesting features of the model and also to show that the Berkeley Madonna implementation presented here is an accurate representation of the Ikeda et al. model, we show that it faithfully reproduces the results of four simulations whose results are shown in the figures of their article. The BM codes used to generate the results of the following simulations are all provided as Electronic Supplementary Material (see Appendix).

Figure 1 shows the results of a simulation of oral water intake (1 l over 5 min) and intravenous infusion of physiological saline; the left panel shows Fig. 10 from the Ikeda article, and the right panel shows results from our BM model, which are clearly a good match to those in their article.

Fig. 1
figure 1

a Simulation of oral water intake (solid lines) and intravenous infusion of physiological saline (dashed lines), both at a rate of 1000 ml per 5 min (see Fig. 10 in Ikeda et al. (1979)). b The same simulations were carried out in Berkeley-Madonna. We simulate, during 3 h, the responses of body fluid and kidney parameters to acute water loading (solid lines) at a rate of 200 ml/min during 5 min (rate of drinking, QIN=0.2 l/min from t = 5 to 10 min) and to intravenous normal saline infusion (dashed lines), solution of 0.9 % w/v of NaCl, containing 154 mEq/l of \(\hbox {Na}^+\) and \(\hbox {Cl}^-\), at the same rate during 5 min (from t = 5 to 10 min, the rate of intravenous water input was QVIN = 0.2 l/min , and intake rate of sodium and chloride was YNIN = YCLI = 30.8 mEq/min). For the simulation of oral water intake (Online Resource 02), the user must replace the following line of BM code: QIN = 0.001 with: QIN = IF (TIME \(\ge\) 5 AND TIME \(\le\) 10) THEN 0.2 ELSE 0.001. For the simulation of intravenous infusion of physiological saline (Online Resource 03), the user must replace the following lines of BM code: QVIN = 0, YCLI = 0.1328 and YNIN = 0.12 with: QVIN = IF (TIME \(\ge\) 5 AND TIME \(\le\) 10) THEN 0.2 ELSE 0, YCLI = IF (TIME \(\ge\) 5 AND TIME \(\le\) 10) THEN 154*0.2 ELSE 0.1328, YNIN = IF (TIME \(\ge\) 5 AND TIME \(\le\) 10) THEN 154*0.2 ELSE 0.12. We observe the rate of urinary output (QWU), the plasma volume (VP), the volume of extracellular fluid (VEC), the intracellular fluid volume (VIC), the plasma osmolality (OSMP), the interstitial fluid volume (VIF), the systemic arterial pressure (PAS), the standard bicarbonate at pH = 7.4 (STBC), the effect of antidiuretic hormone (ADH), and the effect of aldosterone (ALD)

Figure 2 shows the transient response of respiratory parameters to inhalation of 5 % \(\hbox {CO}_2\) over 30 minutes; the left panel shows Fig. 11 from the Ikeda article, and the right panel shows results from our BM model.

Fig. 2
figure 2

a Simulation of the transient response of the respiratory system to 5 % \(\hbox {CO}_2\) inhalation (see Fig. 11 in Ikeda et al. Ikeda et al. (1979)). b The same simulation was carried out in Berkeley-Madonna (Online Resource 04). We simulate, during 1 h, the transient response of the respiratory parameters to the inhalation of 5 % \(\hbox {CO}_2\) in air over 30 min (volume fraction of \(\hbox {CO}_2\) in dry inspired gas FCOI = 0.05 from t = 5 to 35 min). The user must replace the following line of BM code: FCOI = 0 with: FCOI = IF (TIME \(\ge\) 5 AND TIME \(\le\) 35) THEN 0.05 ELSE 0. We observe the alveolar ventilation (VI), the pressure of \(\hbox {CO}_2\) and \(\hbox {O}_2\) in the alveoli (PCOA and PO2A), and the concentration of bicarbonate of the extracellular fluid (XCO3)

Figure 3 shows results from a simulation of glucose tolerance test (infusion of 50 g of glucose over 1 h), including insulin secretion due to a concomitant decrease of extracellular fluid potassium concentration; as above, the left panel shows Fig. 12 from the Ikeda article, and the right panel shows the corresponding results from our BM model.

Fig. 3
figure 3

a Simulation (Fig. 12 in Ikeda et al. Ikeda et al. (1979)) of the glucose tolerance curve with the extracellular fluid potassium concentration. b The same simulation was carried out in Berkeley-Madonna (Online Resource 05). We simulate, during 3 h, a test of glucose metabolism, corresponding to the infusion of glucose at a rate of 1 g/min during 50 min (intake rate of glucose YGLI = 1000 from t = 5 to t = 55 min). The user must replace the following line of the BM code: YGLI = 0 with: YGLI = IF (TIME \(\ge\) 5 AND TIME \(\le\) 55) THEN 1000 ELSE 0. We observe the ECF glucose concentration (XGLE), the ECF potassium concentration (XKE), the plasma osmolality (OSMP), the rate of urinary output (QWU), the renal excretion of glucose (YGLU), and the rate of renal loss of potassium (YKU)

Figure 4 shows, in acid–base disturbances, the central role of the kidney in the compensatory reactions of the body when the normal response of respiration does not occur. The long-term time course of the model behavior in respiratory acidosis or alkalosis is depicted on the pH-[\(\hbox {HCO}_{3}\)] diagram. The response to 10 % \(\hbox {CO}_2\) inhalation and the response to hyperventilation are observed. The right panel shows the results from our BM model, which are in good agreement with the results of Ikeda article, shown on the left panel. The sequence of steps necessary to reproduce this figure with BM implementation is detailed in the specific BM code listing (Online Resources 06 & 07).

Fig. 4
figure 4

a Simulation (Fig. 13 in Ikeda et al. Ikeda et al. (1979)) of respiratory acidosis and alkalosis with renal compensation. Point O shows the normal value of the model of the pH-[\(\hbox {HCO}_3\)] plane. Triangle indicates the plotting of simulated response to 10 % \(\hbox {CO}_2\) inhalation for 48 h, and Filled circle indicates that of hyperventilation, in which VI was fixed at 15 1/min. Equi-pressure lines of \(\hbox {PCO}_2\) are shown with dotted lines for the \(\hbox {PCO}_2\) values of 13.3, 40.0, and 73.0  mmHg. b The same simulations were carried out in Berkeley-Madonna. We first simulate (Online Resource 06), during 48  h, the response to 10 % \(\hbox {CO}_2\) inhalation (volume fraction of \(\hbox {CO}_2\) in dry inspired gas FCOI at the value of 0.1, rather than 0, during the whole simulation and equation (1) unmodified). The bicarbonate concentration of the extracellular fluid (XCO3) and the pH of arterial blood (PHA) are measured at various times from 12 min to 48 h and plotted with Triangle line. We then simulate (Online Resource 07) during 48 h the response to hyperventilation, in which VI was raised to three times normal (alveolar ventilation VI is kept constant to 15  l/min, VI=15, replacing equation (1) of the BM code during the whole simulation). The volume fraction of \(\hbox {CO}_2\) in dry inspired gas FCOI is set at its normal value 0. XCO3 and PHA are measured at various times from 12 min to 48 h and plotted with Filled circle line

4 Discussion

Efforts towards reusability and interoperability have made progress in recent years, not only in the modeling of kidney physiology (Thomas 2009) but also in the wider context of physiology and systems biology (Hunter et al. 2013). For instance, SBML (the Systems Biology Markup language)Footnote 1 (Hucka et al. 2003) is widely used for metabolic networks and models of cell signal transduction, the CellML repositoryFootnote 2 contains several hundred marked-up legacy models (mostly at the level of membrane transport or signal transduction), the JSim Consolidated Model DatabaseFootnote 3 indexes 73390 models across five archives, and annotation tools such as the RICORDOFootnote 4 resource (de Bono et al. 2011) and the ApiNATOMYFootnote 5 (de Bono et al. 2012) project now facilitate the sharing (and even the merging) of physiology and systems biology models.

The present work complements previous re-implementations of the Ikeda model; e.g., a Pascal version was used in teaching at the University of Limburg, Maastricht (Min (1982); Pascal source code in Min (1993)), and extensions of parts of the Ikeda model were used in the Golem simulator (Kofranek et al. 2001). The present Berkeley Madonna version also complements our re-implementations of the early Guyton models (Hernandez et al. 2011; Moss et al. 2012; Thomas et al. 2008) and recent models focused on the kidney itself (Karaaslan et al. 2005, 2014; Moss et al. 2009; Moss and Thomas 2014) or on the role of the kidney in blood pressure regulation (Averina et al. 2012; Beard and Mescam 2012). We provide here a convenient implementation of the Ikeda et al. (1979) model in order to facilitate not only its use in its original form but also to enable its extension. One such improvement would be the incorporation of a more complete model of the RAAS system, which is now much better understood and for which a detailed model has recently been published (Guillaud and Hannaert 2010).