In Silico Biology 6, 0018 (2006); ©2006, Bioinformation Systems e.V.  


A simple mathematical model of adaptation to high osmolarity in yeast


Peter Gennemark1*, Bodil Nordlander2, Stefan Hohmann2 and Dag Wedelin1




1 Department of Computer Science and Engineering, Chalmers University of Technology
SE-412 96 Göteborg, Sweden.
2 Department of Cell and Molecular Biology/Microbiology, Göteborg University
SE-405 30 Göteborg, Sweden.



* Corresponding author

   Email: peterg@chalmers.se





Edited by H. Michael; received December 06, 2005; revised and accepted March 17, 2006; published April 29, 2006



Abstract

We present a simple ordinary differential equation (ODE) model of the adaptive response to an osmotic shock in the yeast Saccharomyces cerevisiae. The model consists of two main components. First, a biophysical model describing how the cell volume and the turgor pressure are affected by varying extra-cellular osmolarity. The second component describes how the cell controls the biophysical system in order to keep turgor pressure, or equivalently volume, constant. This is done by adjusting the glycerol production and the glycerol outflow from the cell. The complete model consists of 4 ODEs, 3 algebraic equations and 10 parameters. The parameters are constrained from various literature sources and estimated from new and previously published absolute time series data on intra-cellular and total glycerol. The qualitative behaviour of the model has been successfully tested on data from other genetically modified strains as well as data for different input signals. Compared to a previous detailed model of osmoregulation, the main strength of our model is its lower complexity, contributing to a better understanding of osmoregulation by focusing on relationships which are obscured in the more detailed model. Besides, the low complexity makes it possible to obtain more reliable parameter estimates.

Keywords: Saccharomyces cerevisiae, osmotic shock, HOG pathway, Fps1, ODEs, model complexity



Introduction

Osmoregulation involves the adaptation of a cell to varying osmotic pressure in order to maintain size and cellular activities [Hohmann, 2002]. The chemical potential of water is central in osmoregulation, and it can be seen as a measure of the effective water concentration in a given area. The water potential is influenced by two factors [Levin et al., 1979]: (1) the osmotic potential and (2) the pressure potential. The first is approximately proportional to the concentration of dissolved molecules of solutes. As the concentration of solute molecules increases, the water potential decreases. The latter takes into account the hydrostatic pressure. If a solution is put under pressure, the water potential increases. For two regions of water with different potentials and separated from each other by a semi-permeable membrane, there will be a water flow to the region of lower potential by osmosis.

Typically, a cell has a higher intra-cellular osmotic pressure (Πi) than extra-cellular osmotic pressure (Πe). The main reason for this difference is that highly charged macromolecules and metabolites attract many small inorganic ions to the cell interior (the Donnan effect, see e. g. Alberts et al., 1994). Due to this difference water will flow into the cell. In isolation, this would cause the cell to swell and potentially lead to cell rupture.

The yeast Saccharomyces cerevisiae prevents the fundamental problem of water inflow and cell swelling by its cell wall, which is less elastic than the plasma membrane. The cell wall resists the expansion of the cell, and creates an inward pressure on the cell contents [Gervais and Beney, 2001]. This pressure is called the turgor pressure (Πt), defined as the difference in the hydrostatic pressure between the inside and the outside of the cell. At equilibrium (equil.), the water potential is equal inside and outside of the cell and the turgor pressure balances the difference in osmotic pressures as [Smith et al., 2000]

Πi = Πe + Πt   (equil.) (1)

A hyper-osmotic shock is a sudden increase in the extra-cellular osmotic pressure, for instance due to the addition of salt to the cell medium. The immediate effect on yeast to an osmotic shock involves water outflow and decreasing volume. The rate of water flow is proportional to the difference between (Πe + Πt) and Πi. In this way, a new equilibrium is reached, in which the higher extra-cellular osmotic pressure is balanced by an increased intra-cellular osmotic pressure (due to the reduced volume), and reduction of turgor pressure (due to reduced size of the cell wall). For a sufficiently large osmotic stress, the turgor pressure is abolished and the volume is significantly reduced at the new equilibrium [Gervais and Beney, 2001; Hohmann, 2002]. We will refer to these processes as the the biophysical system of the cell, see Fig. 1 for an overview.



Figure 1: Overview of osmoregulation in yeast. The biophysical system is illustrated to the right and the control system to the left in the figure. An input stress causes the external osmotic pressure to increase and water to flow out of the cell. In turn, this rapidly leads to turgor pressure and volume decrease. The control system aims at accumulating glycerol in order to regain volume and turgor pressure. The HOG signalling pathway initiates a transcriptional response of several genes involved in glycerol production and the outflow of glycerol is adjusted by Fps1. Arrows indicate dependencies considered in the model.


Generally, the cell strives to keep volume and turgor pressure constant and independent of environmental changes. It therefore has a control system responding to these changes by accumulating glycerol and thereby increasing the intra-cellular osmotic pressure in order to regain volume and turgor pressure [Gervais and Beney 2001; Hohmann, 2002; de Nadal et al., 2002]. This control system consists of two main components as illustrated in Fig. 1. First, the transmembrane glycerol transporter Fps1 closes which promotes the accumulation of glycerol. Second, the High Osmolarity Glycerol (HOG) pathway is activated, which in turn leads to transcription of several genes, among them GPD1 and GPP2, whose proteins are involved in glycerol production. The HOG pathway belongs to the class of Mitogen Activated Protein Kinase (MAPK) pathways that are found in all eukaryotic organisms and are important for transmitting and processing signals from the cell membrane into the cell. Typically, a MAPK pathway consists of a sensing system, a cascade of protein kinases and output systems such as transcriptional regulators. Upon activation, the last kinase in the pathway, Hog1, enters the nucleus and induces transcription.

By combining experimental techniques with mathematical modelling Klipp et al., 2005, proposed a mathematical model of the osmoregulation system. The model includes the HOG pathway, carbohydrate metabolism and glycerol production as well as a biophysical model. In total, the model consists of 35 ordinary differential equations (ODEs) as well as 70 parameters, and it is currently the most detailed mathematical description of the adaptation to high osmolarity in yeast. Such detailed models are important in order to completely understand a particular phenomenon and to ultimately obtain a close match between model and experiment. To a large extent the reaction mechanisms are modelled and, hence, the use of black-box relationships between variables that may have no physical connection is minimised. Besides, building a descriptive model with a high level of detail is important for communicating the system and it also allows a large scope of experiments, such as different gene deletions, to be simulated.

In this article, we present a model that most economically captures the essential physics and biology of osmoregulation. It contains variables for volume, turgor pressure, osmotic pressure and glycerol (and some intermediate variables), and is described by 4 ODEs, 3 algebraic equations and 10 parameters. This model complements the more detailed model of Klipp et al., 2005, and we will refer to the model described here as the simple model and to the model from Klipp et al. as the detailed model.

In general, when compared to a more complex model, a model may be simpler in different ways:

  1. More limited modelling scope. This means that parts or aspects of the model may be absent altogether, or replaced by some black-box description which is not intended to model the internal detail of this function.
  2. Simpler functions and equations, such as for example by assuming linearity. In a model of a metabolic pathway, one could for instance replace all Michaelis-Menten kinetics by bilinear kinetics.

The model we present is simpler with respect to 1 but not 2. This means that to the extent a variable or function is at all considered in the simple model, the model is rich enough to describe it with the same numerical accuracy as in the detailed model. The parameters have been estimated using both new and previously published experimental data on glycerol, and the resulting model appears to explain all known experimental data considered in its modelling scope.

We note that both the simple and the detailed model are dynamic models of osmoregulation, and therefore, simulation of time course experiments can be executed. Naturally, static models could also be considered (as an extreme case of approach 1), and such models would in this sense be even simpler than our model. However, the lack of dynamic information would severely restrict their applicability.

There are several general reasons to construct a basic dynamic model of the osmoregulation system:

Although this list points out some important advantages of simple models even if more complex models are available, we note that different types of models and models with different complexity do not necessarily contradict each other, and are useful as complementary views of the same system.

The rest of this paper is structured as follows. In the next section we present our simple mathematical model of osmoregulation and in section "Parameter estimation" we describe how the parameters of the model were estimated. Furthermore, in section "Simulation results" we compare model simulations with experimental data. We end the paper by a discussion.



Our model

Our simple model consists of two main components. First, a biophysical model describing how the cell is affected by varying extra-cellular osmolarity. We present this model in section "The biophysical model" and note that the same biophysical description is used in the detailed model. The second component is a control model giving a simple description of how the cell controls the biophysical system in order to keep turgor pressure, or equivalently volume, constant. This is done by controlling the glycerol production and the glycerol outflow from the cell. This model is presented in section "The control model".


The biophysical model

We model the biophysical system by considering the dependencies between cell volume (V), turgor pressure (Πt), intra-cellular osmotic pressure (Πi) and extra-cellular osmotic pressure (Πe) by equations derived in thermodynamics. There is evidence that e. g. cell wall thickness, membrane composition and vacuole volume also are affected by osmotic stress, see e. g. Hohmann, 2002. However, the importance of these responses is difficult to judge and these processes are typically non-trivial to include in a model, mainly due to lack of data. Therefore, we disregard them in our modelling. For a more elaborate description of the thermodynamic derivation of the model we refer to Levin et al., 1979.

The flow of water over the cell membrane is proportional to (Πi(t) – Πe(t) – Πt(t)), where t is time, and assuming that the volume (V) only changes due to in- and outflow of water we obtain the following differential equation for V

V'(t) = kp1 ( Πi(t) – Πe(t) – Πt(t) ) (2)

where kp1 is a constant for the hydraulic water permeability constant (usually denoted Lp) times the cell membrane area which we assume constant. At equilibrium, i. e. constant volume and no net flow of water over the membrane, we note that (2) reduces to (1).

We consider glycerol as the sole osmolyte and, hence, ions and other small molecules that have been reported to change upon osmotic shock, see e. g. Sunder et al., 1996, are not considered. This simplification is to a certain extent motivated by experimental results from Reed et al., 1987, who found that glycerol counter-balances in the order of 80% of applied NaCl in S. cerevisiae. For this reason, the intra-cellular osmotic pressure is calculated by van't Hoff's law according to

Πi(t) = ( n + Gly(t) ) / ( V(t) – Vb ) (3)

where Gly [mol] is the main osmolyte glycerol and n are the moles of of other osmotically active compounds in the cell (assumed constant) and Vb is the non-osmotic volume of the cell, corresponding to the sum of the volumes of hydrophobic cellular components (such as lipid bilayers) that are osmotically unresponsive. We see that increased glycerol increases the intra-cellular osmotic pressure. In this way, glycerol can be used to control the turgor pressure of the cell.

The extra-cellular osmotic pressure is only modified by the input signal, for example applied salt stress, and therefore independent of changes in other variables, such as extra-cellular glycerol. This is a reasonable approximation because of two main reasons: first, the intra-cellular volume is much smaller than the extra-cellular volume corresponding to a single cell and second, we do not consider import to the cell (e. g. glucose) and its effects on Πe and it is therefore natural also to leave out the effects of export.

Because of the nature of the elasticity in the cell wall, the turgor pressure is usually dependent on the volume. We consider the approximation of letting turgor pressure be linearly dependent on the volume according to Levin et al., 1979:

Πt(t) = ε ( V(t) / V(0) – 1 ) + Πt(0) (4)

where V(0) and Πt(0) are the initial volume and turgor pressure, respectively, and ε is the so-called volumetric elastic modulus. Restricting Πt to positive values and using the notation VΠt=0 for the volume when Πt = 0, we can rewrite (4) as

Πt(t) = Πt(0) ( V(t) – VΠt=0 ) / ( V(0) – VΠt=0 ) (5)

when V(t) > VΠt=0 and Πt(t) = 0 otherwise.

In the next section we suggest a model for controlling this biophysical model when exposed to extra-cellular osmotic stress.


The control model

To build a model of osmoregulation it is crucial to understand what control systems are available and what type of sensors are employed by these control systems. We consider two parallel ways of control: first, osmotic shock closes the Fps1 channel and thereby allows the cell to accumulate glycerol, and second, osmotic shock triggers the HOG signalling pathway which, in turn, leads to HOG-induced glycerol production after a certain time delay. Concerning Fps1, it is known that the channel closes upon high osmotic pressure, but the sensing mechanism is unknown. For the HOG pathway, there are at least two independent sensors and one of them, Sln1, has been shown to respond to changes in turgor pressure [Reiser et al., 2003], although the exact sensing mechanism is unknown. The other sensor of the HOG pathway (Sho1-branch) is not identified and the mechanism is unknown.

It is believed that it is evolutionary favourable for the cell to strive to keep both turgor pressure and volume reasonably constant. Since it is known that the yeast cell senses turgor pressure by the Sln1 sensor, we use turgor pressure as the controlling variable also in our model.

In order to construct a simple model of how the cell controls the biophysical system (modelled as (2), (3) and (5)), we consider the control problem of keeping the turgor pressure constant at Πt(0) (reference signal), given disturbances in the external osmotic pressure. In practice, both turgor pressure and volume are controlled since volume is linearly dependent on turgor pressure in the biophysical model. See Fig. 2 for an overview of the model.



Figure 2: Overview of the model where the glycerol level is controlled by two parallel ways: (1) adjusted glycerol outflow capacity via Fps1 and (2) increased glycerol production triggered by the HOG pathway.


We define the input signal to the controllers, e, as

e(t) = Πt(0) - Πt(t). (6)

The first control way concerns the adjustment of glycerol inflow/outflow over the Fps1 channel. The Fps1 control function, uFps1, models the glycerol permeability over the channel. The function outputs real values in the interval [0,kp2], where 0 corresponds to completely closed and where kp2 is the glycerol permeability coefficient in a completely open Fps1 channel. For simplicity, we consider proportional control for uFps1 and obtain

uFps1(t) = kp2 ( Πt(0) - e(t) ) / Πt(0) (7)

when e(t) > 0 and uFps1(t) = kp2 otherwise. We model the exchange of internal and external glycerol over the Fps1 channel, uDiff, by Fick's first law of diffusion as

uDiff(t) = uFps1(t) ( Gly(t) / (V(t) - Vb ) - Glye(t) / Ve ) (8)

where Ve is the extra-cellular volume. We consider a single cell in our model and Ve should therefore be interpreted as the fraction of the extra-cellular volume belonging to each cell.

The second way of control concerns the glycerol production induced by the HOG pathway. The control function of the HOG pathway, uHOG, is constrained to positive numbers and assuming proportional control with constant kHOG we obtain

uHOG(t) = kHOG e(t) (9)

when e(t) > 0 and uHOG(t) = 0 otherwise.

We note that we have used very simple proportional controllers, although it is known that for example MAPK signalling pathways often exhibit a switch-like behaviour [Huang and Ferrell, 1996; Kholodenko et al., 1997; Ferrell, 1998]. In principle, we could modify our controllers to some switch-like function, e. g. the logistic function with two parameters, and thereby increase the goodness-of-fit to experimental data. This is natural since the parameter space is enlarged. It is, however, not obvious which switch-like function to select and therefore, in order to avoid an unnecessarily complex model, we use proportional control.

We would also like to point out that phosphatases targeting Hog1 are not explicitly included in the control model, although they are generally believed to control the level of osmotic stress induced phosphorylated Hog1. For instance, there is experimental evidence that phosphatases targeting Hog1 are up-regulated upon osmotic stress [Jacoby et al., 1997; Wurgler-Murphy et al., 1997]. However, given sensors for turgor pressure this up-regulation is superfluous for successful adaptation. Instead, the basal level of phosphatases is enough to dephosphorylate the kinases once the input signal is off [Klipp et al., 2005].

In order to take the time required for transcription and translation into account, we delay the signal from the HOG pathway to glycerol adjustment by a first-order delay approximation using the following ODE

vHOG(t)' = ( uHOG(t) - vHOG(t) ) / td (10)

where vHOG(t) is the time delayed variable and td is the time delay. Alternatively, we could have used an explicit delay with similar results. From a practical point of view, however, (10) is easier to handle since it can be integrated in any standard ODE solver.

The two ways of control are combined in the ODE for intra-cellular glycerol, Gly, as

Gly'(t) = vHOG(t) - uDiff(t). (11)

Similarly, the ODE for extra-cellular glycerol, Glye, is obtained by taking the diffusion term into account as

Glye'(t) = uDiff(t). (12)

Finally, we would like to point out that the use of two parallel ways of control is not unnecessarily complex but rather necessary, in order to adjust the glycerol concentration to a desired value. This is because the flexibility of the controllers are constrained: the HOG controller only assumes non-negative values and the Fps1 controller is dependent on the difference in glycerol concentration over the cell membrane.



Parameter estimation

In order to estimate the parameters of the model we consider both new and previously published data. In particular, we use a data set of absolute time series data for intra-cellular and total glycerol, see Tab. 1. This data set includes seven experiments on wild-type cells with different input signals of NaCl stress. Furthermore, one experiment considers a constitutively open Fps1 mutant.


Table 1: Experiments considered in parameter estimation. The added osmolyte is NaCl in all experiments.
Exp. Input signal Genomic background Var. N Rep. Reference
1 0.17M w.t. G 2 3 Reed et al., 1987
2 0.5M w.t. G
GT
2
5
3
2
Reed et al., 1987
Klipp et al., 2005
3 0.86M w.t. G 2 3 Reed et al., 1987
4 1.0M w.t. GT 5 2 Appendix A
5 0.5M t=0
0.5M t=30
w.t. GT 5 2 Appendix A
6 0.5M t=0
0.5M t=60
w.t. GT 6 2 Appendix A
7 0.5M t=0
0.5M t=160
w.t. GT 7 2 Appendix A
8 0.5M Open Fps1 GT 7 2 Klipp et al., 2005
Abbreviations: Var. = measured variable, N = number of data-points in the time series, Rep. = number of replicates, w.t. = wild-type, G = intra-cellular glycerol concentration and GT = total glycerol concentration. The osmotic pressure (in Osm) of a certain molar of NaCl is obtained by multiplying by 2 (NaCl dissociates into two ions) and then multiplying by the osmotic coefficient of NaCl, 0.93 [Robinson and Stokes, 1959].


The data set of Tab. 1 is incomplete in the sense that only the concentration of glycerol and total glycerol are observed. Besides, data is scarce and monotonic. In order to obtain reliable parameter estimates from such a data set it is useful to consider a simple model and to constrain the parameters. In section "Constraining the model parameters", we present reasonable lower and upper bounds for all parameters. These constraints are obtained not only from the data set of Tab. 1, but also from various literature sources as well as our own data on cell volume immediately after osmotic stress (Appendix B).

We aim at estimating the constrained parameters given the time series data presented in Tab. 1. However, a potential problem in parameter estimation is that all parameters may not be unambiguously determined from the considered data set. Therefore, in section "Observability test", we use a method for algebraic observability to examine whether the parameters of the simple model in theory can be deduced from the experimental data set considered. It turns out that two of the parameters are not observable. However, we overcome this problem by additional reasonable assumptions and, finally, in section "Parameter estimation" we estimate the parameters numerically from the time series data of Tab. 1.

Concerning our data set on intra-cellular glycerol we would like to point out that absolute data from Reed et al., 1987 has been employed, although data with higher sampling frequency is available in Klipp et al., 2005. However, the latter data is relative and the fold increase in intra-cellular glycerol is only in the order of ten for an osmotic shock of 0.5M NaCl. This is not realistic considering that such an osmotic shock should give in the order of 1M intra-cellular glycerol [Oliveira et al., 2003; Reed et al., 1987]. Given only a ten-fold increase in intra-cellular glycerol would indicate an initial concentration in the order of 0.1M, which is about 1000 times greater than values in the literature for non-stressed cells. One possible explanation for these experimental results could be that the cells accidentally have been pre-stressed and, consequently, started to accumulate glycerol before the true stress input signal. A more reasonable explanation is, however, that the method for measuring intra-cellular glycerol is inexact for low concentrations. This is supported by data from Sunder et al., 1996, who use a similar method for measuring intra-cellular glycerol. In their study, the initial concentration varied between 'not detectable' and 0.05 (relative unit), and the concentration one hour after stress was in the order of 0.5 (relative unit).


Constraining the model parameters

The complete model, including the biophysical model and the control model, contains 14 parameters. However, four of these are dependent parameters which we do not need to constrain. Hence, we seek to calculate or constrain the remaining 10 parameters from experimental data. As in the detailed model, we scale the initial volume of the cell to 1. This assures that all variables of the model have values of similar order of magnitude, which facilitates numerical integration. For an overview of all parameters as well as the constraints, we refer to Tab. 2. In this table we also describe how the dependent parameters are calculated.


Table 2: Model parameters and parameter bounds corresponding to lower and upper bounds.
Parameter Bounds
Gly(0) Initial Gly [1.1E-4, 5.0E-4]
Πi(0) Initial Πi [0.60,0.70] Osm-1
Πe(0) Initial Πe [0.24,0.25] Osm-1
Vb Non-osmotic volume [0.31,0.46]
Ve External volume [1.0E3,5.0E3]
kp1 Water perm. coeff. [5.2E-3,160] Osm-1
kp2 Gly perm. coeff. in Fps1 (0,∞)
VΠt=0 V when Πt=0 [0.5,0.99]
kHOG HOG pathway control const. (0,∞) Osm-1
td Time delay [5,30] min
Dependent parameter Value obtained as
V(0) Initial V (relative value) 1
Glye(0) Initial Glye Ve Gly(0) / ( V(0) – Vb )
Πt(0) Initial Πt Πi(0) – Πe(0)
n No. of osmotically active moles other than Gly Πi(0) ( V(0) – Vb ) - Gly(0)
All volumes are scaled such that the initial volume of the cell is 1. Both Gly and Glye represent number of molecules (mol scaled by V(0)).


The parameter constraints were obtained as follows:


Observability test

Using a method for algebraic observability [Sedoglavic, 2002] we examined whether the parameters of the simple model in theory can be deduced from experiments where Πe is input signal and where the concentration of intra-cellular glycerol and total glycerol are output signals. The output signal for glycerol concentration is calculated as

[Gly(t)] = Gly(t) / (V(t) – Vb) (15)

and the corresponding expression for total glycerol is obtained as

[Glytotal(t)] = ( Gly(t) + Glye ) / ( V(t) – Vb + Ve ). (16)

Using these settings, we found that nine of the parameters (all but td) of the model are not observable for these particular in- and outputs. Hence, infinitely many parameter value combinations can fit the observed data. However, by assuming two of the non-identifiable parameters known, we obtain algebraic observability.

For this reason we selected two parameters, Πe(0) and kp1, to be fixed. Concerning Πe(0), we see from Tab. 2 that it to a large extent can be assumed known from its tight bounds. For simplicity, we assign Πe(0) the value calculated from the medium considered in Appendix C. In contrast, for kp1 we cannot rely on tight bounds. However, preliminary simulations of the model with allowed parameter values suggest that the observed data is insensitive to variations in kp1. Therefore, we cannot expect an exact estimate of this parameter to be deduced from our data set. We note that this parameter is also difficult to estimate experimentally, as indicated by the difference of five orders of magnitude in the estimates from the literature. The main reason for this uncertainty is the difficulty of observing the rapid passive response experimentally. Therefore, we assigned kp1 = 1, corresponding to the average value of the lower and upper bounds in logarithmic space. We would like to point out that fixing Πe(0) and kp1 do not significantly affect the estimated values of the other parameters.

It is important to note that a positive algebraic observability test only guarantees that the parameters can be observed when using ideal data. For a realistic data set, which is usually sparse and noisy, additional information might be required for reliable parameter estimates. Concerning our data set for the simple model, we include additional information in form of the parameter bounds as well as the experiment on a modified system (a constitutively open Fps1 mutant).


Parameter estimation

We estimated the remaining eight parameters by minimising the sum of the squares of the difference between simulated and experimental data given by Tab. 1. For one time series experiment the error is calculated as

error = Σi ( X(ti) - Y(ti) )2 (17)

where i indexes the measurement points and where X and Y denote simulated and experimental values, respectively. The total error is obtained by summing over all measured values in all experiments. This error function does not weight the data-points but it is sufficient for our purposes since we only measure glycerol and all data-points have the same order of magnitude. We note, however, that other error functions can be considered. All simulations were run in Matlab (MathWorks Inc.) using the function ode15s for numerical integration.

We searched for a possible global minimum point of the error function by evaluating several randomly chosen starting points in the feasible region. For practical reasons, the upper bounds of infinity were replaced by reasonably large numbers. For parameter sets with sufficiently low error (approximately 0.5% of the trials), the search was continued with the Matlab (MathWorks Inc.) minimisation routine fmincon. The best parameters found are given in Tab. 3. These parameters are obtained also when repeating the parameter estimation, why we are relatively confident that a global minimum of the error function is obtained.


Table 3: Model parameters.
Parameter Value
Gly(0) Initial Gly 2.00E-4
Πi(0) Initial Πi 0.636 Osm
Πe(0) Initial Πe 0.240 Osm
Vb Non-osmotic volume 0.368
Ve External volume 4.79E3
kp1 Water perm. coeff. 1.00 Osm-1
kp2 Gly perm. coeff. in Fps1 0.316
VΠt=0 V when Πt=0 0.990
kHOG HOG pathway control const. 0.416 Osm-1
td Time delay 8.61 min
The values of Πe(0) and kp1 were fixed as described in section "Observability test" while the remaining 8 parameters were estimated from time series data. Gly represents number of molecules (mol scaled by V(0)).


The sensitivity of the error function to perturbations in the parameters was investigated and the result is presented in Tab. 4. The sensitivities indicate reasonable ranges for the parameters. For instance, the parameter kHOG is likely to be found within a factor 0.5 and 2 from its estimated value, because of the steep increase in error for such perturbations. However, we note that this kind of interpretation may not hold if the found minimum of the error function is local and the global minimum is situated in a different region of the parameter space. The sensitivities confirm that variations in kp1 do not have a strong impact on the error function. Actually, the error function is hardly affected for any feasible value of kp1. Similarly, the error function is also insensitive to changes in Vb. However, in this case we have a relatively good knowledge of the parameter value due to the tight bounds. We generally note that parameter identification would have been extremely difficult without the lower and upper bounds. The estimates of the remaining parameters are better, although the estimated parameter values should be interpreted as plausible values.


Table 4: Sensitivity of the error function to variations in the parameters.
Parameter Perturbation
0.1 0.5 0.9 0.99 1.01 1.1 2 10
Gly(0) 1.01 1.000 1.004 1.01 1.81
Πi(0) 1.017 1.025 1.52
Πe(0) 1.002
Vb 1.00 1.000 1.000 1.00
Ve 37.8 1.46 1.004 1.004
kp1 1.00 1.00 1.00 1.000 1.001 1.00 1.00 1.00
kp2 18.9 3.73 1.12 1.006 1.004 1.08 2.15 6.06
VΠt=0 1.38 1.021
kHOG 29.2 8.31 1.31 1.007 1.007 1.23 11.5 120
td 1.00 1.000 1.001 1.01 1.22
The parameters are modified by different perturbations between 0.1 and 10 times the estimated parameter value. Only one parameter is modified at a time. The values in the table correspond to the error of the perturbed model divided by the error of the model with estimated parameters. When a parameter value is outside its bounds the model error is set to infinity.



Simulation results

In this section we present simulated experimental data using the simple model. First, in section "Wild-type strain with NaCl input stress" we compare simulated data with some of the time series experiments of Tab. 1. Because of the relatively few available time series, we used all these data for parameter estimation and no such data remains for model validation. Instead, we validate the model by qualitative comparisons with mutant experiments considered in Klipp et al., 2005, and with mutant experiments with different polyol input signals studied by Karlgren et al., 2004. These two comparisons are presented in section "Mutant strains with NaCl input stress" and section "Mutant strains with various polyol input stress", respectively.


Wild-type strain with NaCl input stress

Simulated and experimental data for some of the wild-type experiments of Tab. 1 are given in Fig. 3. We note that the model parameters have been estimated using all experiments, and hence, the goodness-of-fit for an individual experiment is not optimised.



Figure 3: Simulation of osmotic shock of NaCl in wild-type cells. Upper plot: 0.5M NaCl, middle plot: two osmotic shocks (t = 0 and t = 30) of 0.5M NaCl, lower plot: two osmotic shocks (t = 0 and t = 60) of 0.5M NaCl Experimental data for total glycerol is indicated by squares and data for intra-cellular glycerol is indicated by triangles. In the upper plot, time-series of intra-cellular and total glycerol from a corresponding simulation of the detailed are also depicted.


The upper plot of Fig. 3 depicts experiment 2 of Tab. 1. The input to this experiment is 0.5M NaCl, corresponding to an increase in the extra-cellular osmotic pressure by 0.93 Osm. Experimentally, glycerol and total glycerol have been measured and these data are plotted together with simulated data from the model. Simulated data illustrates how volume and turgor pressure drop immediately upon osmotic stress and subsequently regain due to glycerol accumulation. Turgor pressure, which is the controlled variable, is not increasing to its original value during the time span of the experiment. The volume, on the other hand, is almost completely recovered to its initial value after a couple of minutes. The main reason for incomplete recovery is that the model parameters are estimated using data, where the measured glycerol concentration is not sufficient for complete recovery of both volume and turgor pressure. The reason why volume and not turgor pressure is recovered is the estimated value of the parameter VΠt=0, indicating a low elasticity of the cell wall, and no turgor recovery until volume is almost completely recovered. For a lower value of VΠt=0 (higher elasticity), the turgor pressure would recover faster while the volume would recover slower.

An interesting aspect is how well the simple model behaves in comparison with the detailed model. However, a quantitative direct comparison between the two models is not easily done without additional data. The models have been fitted with different (only to a limited extent over-lapping) data sets, and it would not be not fair to judge the quality of one model by predicting the data that is used to fit the other model. In an attempt to at least illustrate some conceptual differences between the models, we have included simulations of intra-cellular and total glycerol from the detailed model in the upper plot of Fig. 3. Notably, the detailed model fails to accurately capture the levels of total glycerol. For further discussion of this we refer to section "Discussion".

Experiments 5 and 6 of Tab. 1 are given in the two lower plots of Fig. 3. Here, two subsequent osmotic shocks of NaCl are applied and simulated data is compared to experimental data of total glycerol. These simulations suggest that the model is able to control series of osmotic shocks, something that was experimentally demonstrated in Klipp et al., 2005.


Mutant strains with NaCl input stress

The model includes two parallel ways of control: the glycerol outflow control via Fps1 and the glycerol production initiated by the HOG signalling pathway. An interesting aspect studied by Klipp et al., 2005, is how the system and the model behave with only one way of control.

The importance of the Fps1 control function can be studied by a mutant experiment with constitutively open Fps1, an experiment that is also included in our set of training experiments. To simulate this experiment, we manually fix uFps1 = kp2 and we also adjust Gly(0) to obtain a realistic initial value of total glycerol. Experimental data suggests an over-production of glycerol and a prolonged activation time of the HOG pathway. As illustrated in the upper part of Fig. 4, simulated data from the simple model is in accordance with these findings. In particular, the glycerol production is almost doubled compared to the wild-type experiment. The prolonged activation time of HOG signalling cannot be explicitly observed since Hog1 is not a variable in the simple model. However, the time span of turgor pressure equal to zero gives an implicit measure of the time span of Hog1 activation.



Figure 4: Simulation of one osmotic shock of 0.5M NaCl in modified systems. Upper plot: a mutant with constitutively open Fps1, middle plot: a gpd1Δgpd2Δ mutant, lower plot: a mutant with constitutively open Fps1 and over-expression of phosphatases. Experimental data for total glycerol is indicated by squares.


We can also test the importance of the second way of control by considering a gpd1Δgpd2Δ double mutant unable to produce glycerol. To simulate this experiment, we assign kHOG = 0 and set the initial concentrations of glycerol to zero. See the middle part of Fig. 4 for simulated time series data for some of the model variables. Since we have no control of glycerol production, turgor pressure and volume remain low for the complete simulation. Also this is in accordance with experimental observations [Klipp et al., 2005].

In the same way, we can test the other mutants considered in Klipp et al., 2005. For instance, we can simulate an experiment of over-expression of phosphatases targeting signalling molecules in the HOG pathway in combination with constitutively open Fps1. Experimentally this leads to even slower adaptation than for the experiment with a constitutively open Fps1 and normal phosphatase activity. This is also confirmed in our simulations, see the lower part of Fig. 4. Over-expression of phosphatases was implemented by reducing kHOG by a factor of 2.

These results of mutant strains with NaCl input stress confirm that two ways of control exist in the cell: one increases the concentration of glycerol and the other adjusts the outflow of glycerol. In this way, the simple model strengthen the hypothesis in Klipp et al., 2005, that at least two control functions are necessary to efficiently counter-balance an osmotic shock in the cell.


Mutant strains with various polyol input stress

Karlgren et al., 2004 investigate the osmoregulation system in a gpd1Δgpd2Δ double mutant with constitutively open Fps1 and with input stresses of different polyols, including glycerol. The experimental result with glycerol as input signal suggests a rapid recovery of turgor pressure and volume after stress. This is due to diffusion of glycerol over Fps1. In this way intra-cellular osmotic pressure is immediately balanced and, consequently, the cell does not need to initiate increased glycerol production. This result is in agreement with our simulation of the same system, see the upper plot of Figure 5.



Figure 5: Simulation of one osmotic shock of different stress agents in a gpd1Δgpd2Δ mutant with constitutively open Fps1. Upper plot: 2M glycerol, middle plot: 1M xylitol and lower plot: 1M sorbitol. Note the different time scales between the upper and the two lower plots.


Similarly, Karlgren et al., 2004, stress the cell with xylitol or sorbitol that are both known to diffuse over the cell membrane, but with a slower rate than glycerol. Experimental data suggests that the level of HOG activation is dependent on the diffusion rate: for a slowly diffusing osmolyte like sorbitol, the activation is more pronounced than for xylitol or glycerol. To simulate these experiments we assume that xylitol and sorbitol can diffuse over the open Fps1 and we decrease kp2 by a realistic factor according to experimental data on initial uptake rates in Karlgren et al.. For xylitol we use a factor of 5 for sorbitol a factor of 15. Fig. 5 illustrates that also these simulations are in agreement with experimental results in Karlgren et al., 2004. In particular, the increasing phosphorylation of Hog1 (indicated as the time span of turgor pressure equal to zero) is dependent on the diffusion rate of the considered osmolyte.



Discussion

We have presented a simple ODE model of the adaptation to high osmolarity in S. cerevisiae. The model includes two parallel ways of cellular control: one way controlling the production of glycerol and the other controlling the outflow of glycerol. Our simple model is based on the same basic assumptions as the detailed model [Klipp et al., 2005]. The parameters of the model were obtained by first constraining the values from various data, and then numerically estimating the parameters with respect to absolute time series data for wild-type cells as well as for one genetically modified strain. The complete parameter estimation process is described in detail.

We have validated the model by predicting the behaviour of other genetically modified strains and input functions. In particular, the simple model qualitatively predicts the mutant experiments in Klipp et al., 2005, as well the experiments with different input functions in Karlgren et al., 2004. Therefore, our results indicate that we have constructed a simple but functional model of osmoregulation in yeast.


Comparing the two models

We have chosen to model the essential parts of the system. Our general reasons for this were stated in the introduction and here we discuss these more thoroughly with respect to our model (see also Brooks and Tobias, 1996).

In the simple model, the two parallel ways of control of the osmoregulation system can easily be followed (see e. g. Figure 2), while they are much harder to trace in the detailed model (see e. g. Figure 1 in Klipp et al., 2005). In addition, this simplicity makes it easier to construct and validate the model. What the simple model is able to do with respect to its modelling scope, it does with similar precision as the detailed model.

While the simple model gives a more global view of the system, a consequence is that mechanisms do not have a direct biological counterpart to the same extent as in the detailed model. For instance, the HOG controller in the simple model does not have a direct correspondence in the real biological system. Instead, in combination with the time delay, the HOG controller corresponds to both the HOG signalling pathway, transcription/translation and the synthesis of enzymes involved in glycerol production. Therefore, the HOG control signal can only be qualitatively compared to experimental data. Notably, by extending the model by variables that are accessible for measurement, like Hog1 and Hog1-induced mRNA, we could do these kind of comparisons at the cost of increased model complexity.

The simple model also has a potential for more accurate simulations and predictions, since it is easier to reliably fit the model to currently available data. Roughly speaking, we have more data per parameter, which is natural since we have focused on the most essential parts of the system. The difficulty of the detailed model in predicting an experiment is illustrated in Figure 3, and confirms the difficulty of accurately fitting the detailed model with limited data. Since data in systems biology is often very limited compared to the complexity of the systems, this is a fundamental issue.

Another general observation is that a simple model is usually a good starting point when building more complex models. It is then natural to incrementally add details to the model and simultaneously matching the complexity of the model with both the purpose of the model and available data. In our case, the simple and the detailed models have been constructed in parallel. Since the main characteristics of the system could more easily be observed in the simple model and since the simple model could be parameterised with higher confidence than the detailed model, insights gained from the simple model were useful in constructing the detailed model. For instance, observations on the simple model suggested how to adjust the detailed model to give realistic output on intra-cellular glycerol. Therefore, we believe that the simple model can continue to be a valuable tool for testing new experimental scenarios as well as to constrain further developments of the detailed model.

An argument in favour of more complex models is that the scope of potential system modifications is greater than for simpler models. This is because it is easier to exactly map a certain biological modification into the equations and/or parameters of a complex model. For instance, it is relatively evident how to implement over-expression of phosphatases in the detailed model where phosphatases are explicitly modelled, while it is not straightforward in the simple model where phosphatases are only implicitly modelled.

We note that there are explicit measures of model complexity [Myung and Pitt, 2004; Crampin et al., 2004]. In principle, a measure of model complexity could be used to seek the simplest possible model capturing the trends in data. For instance, we could ask whether we would obtain a better model by removing the time delay from the simple model.


An engineering control model

We will finally discuss how the biophysical model in principle could be controlled using engineering means. Here we wish to point out that these means are not available in the cell and that this scenario is purely artificial. In contrast to the simple model with two parallel ways of control, we now employ a single controller, but let that controller be flexible, in the sense that the glycerol level can be modified to any desired value. As for the HOG controller in the simple model, we let the control signal be time delayed. We employ a PID-controller (Proportional - Integral - Derivative), which is a standard controller in engineering applications. We refer to Fig. 6 (upper) for an overview of the model.

To test the model, we first only consider P-control and manually select a relatively low control coefficient in order to avoid overshoot and oscillations. This leads to a slow adaptation, see the upper plot of Fig. 6. We note that this model is simpler than the simple model, but has less explanatory power since it does not describe the two parallel ways of control.



Figure 6: Controlling the biophysical model by engineering means. The glycerol level is controlled by a standard PID-controller. The control signal, u(t) is derived as u(t) = Kp e(t) + Kie(t) dt + Kd e'(t) where Kp, Ki and Kd are parameters and e(t) is defined in (6). For the time delayed glycerol production we consider Gly'(t) = u( t - td ) and for the biophysical model we use (2), (3) and (5). Parameters are taken from Tab. 3. Upper plot: Simulation of one osmotic shock with Kp = 0.035 and Ki = Kd = 0. Lower plot: Simulation of one osmotic shock with Kp = 0.055, Ki = 0.0001 and Kd = 0.3. For comparison, data for turgor pressure of the simple model is included in both plots.


In order to obtain faster adaptation we consider PID-control and manually select control coefficients such that turgor pressure is recovered rapidly with little oscillation or overshoot, see the lower plot of Fig. 6. This model has a higher complexity than the P-control model, but does not add any explanatory power to the model, since one controller is still employed.

Compared to the simple model we reach reference turgor pressure much faster in the engineering model. The main reason why turgor pressure is only slowly recovering in the simple (or detailed) model is that glycerol diffuses via Fps1 to the medium. On the other hand, turgor pressure is increasing from zero level earlier in the simple model than in the engineering model. One reason for this difference is that we use a delay approximation in the simple model while the engineering model is explicitly delayed. The main reason, however, is that the control coefficient of the simple model (kHOG) is about ten times higher than the proportional control coefficient of the engineering model. The simple model can afford having a high control coefficient since the risk of overshoot and oscillations is largely reduced by the outflow of glycerol via Fps1. In contrast, for the engineering model such a high control coefficient would lead to an overshoot and subsequent oscillations.

Finally, from this engineering approach we would like to emphasise that control of the biophysical system can be obtained given an ideal turgor sensor and perfect but time delayed control of the glycerol level. Although extremely simplified, we believe that this engineering model can be used as a first test for new scenarios that one wants to consider experimentally. One example could be experiments with different osmotic stress input signals.



Authors' contribution

PG conceived the modelling concept, constructed the model, retrieved data, did implementations and analyses and drafted the article.
BN designed, executed and interpreted the experiments and contributed biological background information.
SH supervised the experimental study and the interpretation of biological information.
DW critically edited the manuscript content.



Acknowledgements

Thanks to Gunnar Cedersund (Fraunhofer-Chalmers Research Centre) for suggesting the observability test and to Edda Klipp (Max-Planck Institute for Molecular Genetics, Berlin) for comments on the model.

This work was supported by Chalmers Bioscience Program (PG and BN were partners of a systems biology pair-student project) and by the European Commission via contract LSHG-CT2003-530203 (the QUASI project, to SH).




References



Appendix


A. Experimental data for total glycerol

Total glycerol was measured as described in Klipp et al., 2005.

Table A1: Total glycerol measured for different input signals.
Experiment Total glycerol (mM)
  t=0 min 30 60 90 120 150 180 min
1M NaCl 0.43 0.69 1.21 2.35 2.76    
  0.50 0.63 1.08 1.81 2.60    
Double stress, 30 min 0.46 1.08 1.89 2.39 2.82    
  0.50 1.16 1.86 2.23 2.91    
Double stress, 60 min. 0.35 1.11 1.71 2.43 2.95 3.68  
  0.55 1.18 1.77 2.37 2.98 3.29  
Double stress, 160 min. 0.38 1.13 1.86 1.92 2.60 3.01 4.52
  0.41 0.87 1.69 2.19 2.43 2.46 3.21
Double stress means addition of 0.5M NaCl at t=0 and a second stress of 0.5M NaCl at indicated time points. Each experiment is performed in two repetitions.


B. Relative cell volume immediately after an osmotic shock

To estimate the relative cell volume immediately after an osmotic shock we measured the cell area in ordinary microscope (Leica FW 4000) images and assumed spherical cells. We used the same yeast strain and growth conditions as in Klipp et al., 2005.

Table A2: Relative cell volume immediately after an osmotic shock.
Applied NaCl stress Measured cells Volume
Molar Osm   Mean St.dev.
- - 384 1.00 0.36
0.5 0.93 420 0.70 0.23
1.0 1.86 365 0.61 0.22
1.5 2.79 178 0.53 0.17
The relatively large standard deviations reflect the natural cell size distribution.


C. Estimation of Πe(0)

20 g/l C6H12O6 → 20/180 M = 0.111 M → [i = no of moieties into which the solute dissociates=1, Φ = osmotic coeff. = 1.01] → Π = 0.112 Osm.

5 g/l (NH4)2SO4 → 5/132 M = 0.0379 M → [i =3, Φ = 0.767] → Π = 0.087 Osm.

3 g/l KH3PO4 → 3/136 M = 0.0220 M → [i = 2, Φ = 0.901] → Π = 0.040 Osm.

0.5 g/l MgSO4*7H2O → 0.5/246 = 0.00203 M → [i = 2, Φ = 0.606] → Π = 0.0025 Osm.

Amino acids and vitamins are neglected. Values for the osmotic coefficients are taken from Robinson et al., 1959. The total osmolarity is approximated to 0.24 Osm by summing the individual contributions.