Skip to main content

Deterministic mathematical models of the cAMP pathway in Saccharomyces cerevisiae

Abstract

Background

Cyclic adenosine monophosphate (cAMP) has a key signaling role in all eukaryotic organisms. In Saccharomyces cerevisiae, it is the second messenger in the Ras/PKA pathway which regulates nutrient sensing, stress responses, growth, cell cycle progression, morphogenesis, and cell wall biosynthesis. A stochastic model of the pathway has been reported.

Results

We have created deterministic mathematical models of the PKA module of the pathway, as well as the complete cAMP pathway. First, a simplified conceptual model was created which reproduced the dynamics of changes in cAMP levels in response to glucose addition in wild-type as well as cAMP phosphodiesterase deletion mutants. This model was used to investigate the role of the regulatory Krh proteins that had not been included previously. The Krh-containing conceptual model reproduced very well the experimental evidence supporting the role of Krh as a direct inhibitor of PKA. These results were used to develop the Complete cAMP Model. Upon simulation it illustrated several important features of the yeast cAMP pathway: Pde1p is more important than is Pde2p for controlling the cAMP levels following glucose pulses; the proportion of active PKA is not directly proportional to the cAMP level, allowing PKA to exert negative feedback; negative feedback mechanisms include activating Pde1p and deactivating Ras2 via phosphorylation of Cdc25. The Complete cAMP model is easier to simulate, and although significantly simpler than the existing stochastic one, it recreates cAMP levels and patterns of changes in cAMP levels observed experimentally in vivo in response to glucose addition in wild-type as well as representative mutant strains such as pde1Δ, pde2Δ, cyr1Δ, and others. The complete model is made available in SBML format.

Conclusion

We suggest that the lower number of reactions and parameters makes these models suitable for integrating them with models of metabolism or of the cell cycle in S. cerevisiae. Similar models could be also useful for studies in the human pathogen Candida albicans as well as other less well-characterized fungal species.

Background

Cyclic adenosine monophosphate (cAMP) is an important signalling and regulatory molecule. In eukaryotes cAMP activates Protein Kinase A (PKA), the target kinase of the cAMP-mediated signal transduction pathway. In the widely used model baker's yeast Saccharomyces cerevisiae, this pathway regulates a variety of cellular processes, including metabolism [1], response to stress [2, 3] and progression through the cell cycle [4, 5]. The pathway is modulated by external nutrients, most notably glucose [6]. The transition to growth on glucose in yeast is orchestrated by a tightly regulated pattern of changes in cAMP levels as a result of series of interactions involving the components of the cAMP/PKA pathway (Figure 1). Cyclic AMP is synthesized by adenylate cyclase (Cyr1p), which in turn is regulated by Gpa2p [7] and Ras2p [8], both of which are G proteins. Gpa2p is activated by the G-protein-coupled receptor Gpr1p, which in turn is activated by glucose [9]. Gpa2p is deactivated by the regulator of G protein signalling protein (RGS) Rgs2p, as well as its own intrinsic GTPase activity [10]. Ras2p is activated by the guanine-nucleotide-exchange factor (GEF) Cdc25p [11] and Sdc25p [12], and deactivated by the GTPase activating proteins (GAPs) Ira1p and Ira2p [13]. The level of intracellular GTP is believed to influence the level of GTP-bound Ras2p [14], and the GTP level increases following a pulse of glucose [13], although the mechanism behind this increase is not fully understood.

Figure 1
figure 1

Schematic representation of elements of the cAMP pathway in S. cerevisiae.

The cAMP level is modulated by two phosphodiesterases: Pde2p has higher affinity for cAMP (around 1 × 10-3 mM) [15] compared to Pde1p which has a lower affinity for cAMP in crude extracts (around 0.1 mM) [16, 17]. Yeast cells previously starved for glucose exhibit a characteristic "spike" of cAMP following addition of glucose to the growth media. In wild-type cells, this spike reaches a peak at around 60 seconds before reaching a steady level after around 120 seconds

In the yeast cell, the only known function of cAMP is to activate protein kinase A (PKA). A molecule of PKA consists of two regulatory (R) and two catalytic (C) subunits. Under low cAMP concentrations, the R and C subunits are bound together to form a catalytically inactive heterotetramer. The complex is activated when two molecules of cAMP bind to each R subunit, causing their dissociation from the catalytic subunits. Following dissociation, the free C subunits can phosphorylate their targets. In yeast, the R subunit is encoded by BCY1, while the C subunits are encoded by the partially redundant genes TPK1, TPK2 and TPK3. Recently specific as well as common phosphorylation targets of the Tpk isoforms have been identified [18].

PKA exerts feedback on the system in several ways. First, it has been shown that the low affinity cAMP phosphodiesterase Pde1p is phosphorylated following a glucose pulse and Pde1p can be phosphorylated by bovine PKA [19]. Phosphorylation of Pde1p leads to increased phosphodiesterase activity, which plays a part in reducing the cAMP level following a glucose induced spike. Secondly, PKA can phosphorylate Cdc25p, leading to its dissociation from Ras2p [20]. This results in a decrease in adenylate cyclase activity. Finally, PKA may be able to regulate itself, as it has been demonstrated that Tpk1p is phosphorylated following a glucose pulse [21].

The roles of certain components of the cAMP pathway are still disputed. One of them is that of the Kelch Repeat Homologue proteins Krh1 and Krh2, also known as Gpb1 and Gpb2, as they are believed to function as beta subunits of Gpa2p. According to Harashima and Heitman [22] the Krh proteins stabilize the Ira proteins, the GTPases of the Ras proteins. Deletion of the Krh proteins leads to a loss of the Ira proteins, and therefore cAMP signalling is increased. However, there is evidence that shows that the Krh proteins enhance the association between the regulatory and catalytic subunits of PKA, and this enhancement is removed when the Krh proteins form a complex with activated Gpa2 [7, 23]. Further evidence for the role of the Krh proteins comes from studies of adenylate cyclase (cyr1Δ) mutants [7]. Yeast cyr1Δ pde2Δ mutants can survive on YPD supplemented with 5 mM cAMP. However, the quadruple cyr1Δpde2Δkrh1Δkrh2Δ mutants survive in the presence of 1 mM cAMP, suggesting that the Krh proteins directly inhibit PKA activity, as PKA activity is necessary for yeast survival. In addition a cyr1Δpde2ΔGPA2Q 300Lmutant (with Gpa2 locked in its constitutively active GTP bound state) requires 1 mM cAMP for survival. This gives further support to the theory that Krh is recruited to active Gpa2.

The reductionist approach [24] has taught us much about individual elements of the cAMP pathway; however a quantitative and integrated mathematical representation is needed to fully understand its dynamics. Models of two broad categories can be used for this purpose: deterministic and stochastic [25–27]. Deterministic models which usually consist of a series of ordinary differential equations (ODEs) to describe the system in respect to time, have been used to study yeast systems such as glycolysis [28], the pheromone pathway [29–31] and the cell cycle [32]. Stochastic models on the other hand are used when intrinsic noise is important to the system, such as when low species numbers are involved [33]. However, stochastic models can be computationally expensive to simulate [34].

A stochastic model has been developed to examine the effects of altering the intracellular GTP levels on the Ras/cAMP/PKA pathway [14]. However, in yeast the components of the cAMP pathway are present in high numbers (proteins in thousands, nucleotides in millions) making a deterministic model more appropriate. Moreover, this stochastic model did not include the Krh proteins. In this study we present a deterministic mathematical model of the yeast Ras/PKA/cAMP pathway, with components such as the Krh proteins that have not been included before. Our model has been fitted to experimental data. It is much easier to simulate than is the previously reported stochastic model, yet it can faithfully replicate intracellular species concentrations observed at steady state, and following a perturbation of the system with glucose.

Methods

ODE models of biochemical systems consist of variables and parameters. The variables represent species concentrations, whereas the parameters include rate coefficients, kinetic parameters, etc. If we represent the variables (x i ) as a vector X:

(1)

and the parameters (k i ) as a vector θ:

(2)

then an ODE model can be represented with the following equation:

(3)

The models generated in this study are summarised in Table 1. The reaction formulae which form the basis of the models were entered into Gepasi [35] and/or Copasi [25], and these programs were used for earlier inspection of the models. The models were later exported in Systems Biology Markup Language (SBML) format [36], which allowed the models to be exchanged between programs. SBToolbox in Matlab [37] was used for parameter estimation, parameter sensitivity analysis and model simulations.

Table 1 Summary of the models generated in this study.

Steady state parameter sensitivity analysis was carried out according to the following equation:

(4)

where S ij is the sensitivity of species i in relation to parameter j, Xss i is the steady state level of species i, p j is the value of parameter j, and Δp j is a perturbation of parameter j (equal to 1% of the parameter value).

Cyclic AMP time course data were taken from the literature [19, 38]. As cAMP levels are often reported in terms of nanomoles per gram of wet weight (or equivalent) it was necessary to convert them to nanomolar using the following formula:

(5)

where C(nM) is the nanomolar concentration of cAMP, C(nmolesgww-1) is the cAMP concentration in nanomoles per gram of wet weight reported in the literature, Cw is the conversion factor from grams wet weight to grams dry weight (0.15) and Vc is the volume of 1 × 107 cells in litres (2.68 × 10-6, there are approximately 1 × 107 cells in 1 gram of dry weight).

We recognise that ODE models of this type assume that all cells are identical, which may well not be the case [39].

Parameter estimation

The values of system parameters which were not experimentally derived, were fitted to experimental cAMP time course data using simulated annealing [40, 41], an estimation method that is very efficient in finding a close approximation of the global minimum of an optimization problem. It is based on a probabilistic search, in which every iteration of the algorithm replaces the current solution by a random nearby solution, using a probability distribution that tends to move the solution towards the global minimum. The simulated annealing algorithms found in SBToolbox in Matlab [37] with the SBToolbox function SBparameterEstimation were used for parameter estimation in the current study.

Results

The Protein Kinase A module

The only known biochemical role of cAMP is to activate PKA. This process has a complicated reaction scheme, which is challenging to model. A general guiding principle when building models is to make the model as simple as possible, while capturing realistic behaviour [42]. The expected behaviour of any PKA model must be consistent with the currently available experimental evidence. Firstly, a degree of PKA activity is required for cell viability [43]. If no cAMP is present, the cell is nonviable [44]; therefore all catalytic subunits must be contained within the inactive tetramer in the absence of cAMP. The level of free catalytic subunits must be sensitive to the level of cAMP. The cAMP level can range from 0.015 mM in glucose starved cells, to approximately 0.05 mM (a peak of cAMP induced by a glucose pulse) [38].

The stochastic PKA module reported by Cazzaniga et al [14] makes several assumptions. The binding constants for the association of a cAMP molecule with the PKA tetramer are the same for all cAMP bound states of PKA, as well as the dissociation constants. The underlying assumption is that cAMP binds to PKA in a non-cooperative manner, i.e. the binding of a molecule of cAMP to PKA does not affect the binding/dissociation of further cAMP molecules. In addition, the dissociation of the cAMP-bound PKA holoenzyme, and the subsequent dissociation of cAMP from the free R subunit is considered to be very fast, as is the reassociation of the PKA holoenzyme. We have adopted the same assumptions for our deterministic model.

The stochastic PKA module found in [14] can be converted into a series of deterministic ODEs to give PKA Model A (Table 1). The reactions of this module are summarized in Table 2. The kinetic rate constants are taken from the stochastic time constants found in [14]. This deterministic model can be tested by simulating it over a 100-second time course. Initially the cAMP level is set to 0 and the PKA level is set to 2500 molecules per cell. After 10 s the cAMP level is set to 270900 molecules per cell (equivalent to 0.015 mM). After 30 seconds the cAMP level is increased to 903000 (equivalent to 0.05 mM), and after 60 seconds the cAMP level is decreased to 270900 molecules per cell. For the cAMP level to affect the greatest control on the system, the difference in the level of free catalytic subunits of PKA between low and high cAMP levels should be as high as possible. Cyclic AMP activates PKA, therefore we expect to see an increased difference between active and inactive PKA when cAMP levels are physiologically high.

Table 2 Reactions of PKA Model A

As shown by the blue trace of Figure 2 (panel A) no free catalytic subunit is present when cAMP is set to zero. The model shows changes in the proportion of free catalytic subunits of PKA when cAMP is set to low (C low ) and high (C high ) levels. However the difference between the two states is not great – 27.7% when cAMP is low compared to 40.6% when cAMP is high. It is therefore important to optimize the model, and for this purpose parameter sensitivity analysis was carried out. As shown in Figure 2 (panel B), the parameter k cAMPgain is the most sensitive to variations in PKA level. The parameters of this model were scanned further to identify those which determined the highest difference between C low and C high . Figure 2 (panel C) shows how the difference between C low and C high depends on the parameters k cAMPgain and k cAMPloss. The peak values of this distribution were used to create an optimised model, named PKA Model B, whose simulation is shown by the red trace of Figure 2 (panel A). In PKA Model B, the level of C low now stands at ~10% whilst that of C high is approximately 90%.

Figure 2
figure 2

Deterministic model of the PKA module. (A) Simulation of PKA Model A (blue trace) and PKA Model B (red trace). The cAMP level is 0 initially, and is increased to 270900 molecules per cell (equivalent to 0.015 mM) after 10 seconds, increased to 909000 molecules per cell after 30 seconds, and decreased to 270900 molecules per cell after 60 seconds. (B) Steady state parameter sensitivity analysis carried out on the PKA module. (C) Parameter scan of PKA Model A. The greatest value for PKA difference (79.1%) is achieved when k cAMPgain = 0.1, k cAMPloss = 2.2 × 105, k PKAdiss = 1 × 105, k RcAMPdiss = 100, k PKAass = 1000.

With regards to modelling PKA activation, we wanted to test if the multi-reaction module used so far could be approximated with mass action kinetics or Michaelis-Menten type kinetics. We therefore generated two new PKA Models named C and D, respectively for mass action and Michaelis-Menten kinetics. In these new models, the reaction scheme becomes:

These models are defined by the following ODEs.

PKA Model C:

(6)

PKA Model D:

(7)

We found that these simplified PKA modules could accurately approximate species levels of the optimized PKA Model B, with the following parameters: for PKA Model C, kA = 8.72e-17; kR = 1000; for PKA Model D, V max f = 1e-13; K M F = 1e7; V max r = 1000; K M r = 0.01 (Figure 3).

Figure 3
figure 3

Optimisation of the PKA model. The blue trace shows the simulation of PKA Model B, the red trace – PKA Model C, and the green trace – PKA Model D. The cAMP level is 0 initially, and is increased to 270900 molecules per cell (equivalent to 0.015 mM) after 10 seconds, increased to 909000 molecules per cell after 30 seconds, and decreased to 270900 molecules per cell after 60 seconds.

We also compared steady state proportions of free catalytic subunit of PKA (C free ) of each PKA model as a function of the cAMP concentration (Figure 4). At low cAMP concentrations, the Michaelis-Menten based model (PKA Model D) slightly over-estimated, while the mass action based model (PKA Model C) slightly underestimated the level of C free , respectively, in comparison to the optimised PKA Model B.

Figure 4
figure 4

Steady state levels of free C in the PKA models under various cAMP levels. The blue trace shows the simulation of PKA Model B, the red trace – PKA Model C, and the green trace – PKA Model D. Parameters of PKA Model B are the same as in Figure 3. Parameters of PKA Model C are: k A = 8.72 × 10-17; k R = 1000. Parameters of PKA Model D are: k cat = 10-13; K mF = 107; V maxR = 1000; K mR = 0.01.

The results of simulating these models show that it is possible to simplify the PKA module greatly without loss of performance. It is preferable to use the mass action based module, as it has just three state variables and two parameters. This compares favourably to the complex PKA module which has nine state variables and four parameters. Therefore we adopted the mass action based module to construct the model of the entire cAMP pathway.

Development and simulation of a conceptual model of the complete cAMP pathway

As a step towards developing a deterministic model of the complete cAMP pathway, we first constructed a conceptual model named Simplified cAMP Model A (Table 1). It consists of three ODEs (equations 8–10), based on mass action kinetics and uses unitless species concentrations and parameter values shown in Table 3. Equation 8 represents the combined G proteins activation and inactivation module, equation 9 – the PKA module, and equation 10 – cAMP synthesis and degradation.

Table 3 Parameters of Simplified cAMP Model A
(8)
(9)
(10)

The conservation relationships in Simplified cAMP Model A are described below in equations 11–12 and represent the conservation of the total number of G proteins (GP), and that of total PKA molecules, respectively:

(11)
(12)

where GPi and GPa are the numbers of inactive and active G proteins, respectively. PKAi and PKAa are the number of inactive and active PKA molecules, respectively.

The Simplified cAMP Model A was tested to see if it could reproduce the changes in cAMP levels observed experimentally during a glucose pulse. In a preliminary step, the initial concentrations of cAMP, PKAi and GPi were set to 1, and PKAa, GPa and glucose were set to zero. A steady state was then found, and subsequently all concentrations were set to their steady-state level. A simulation of the model was then run changing the glucose concentration to 5 after 5 time units. As shown in Figure 5 (panel A), a spike of cAMP was observed when the glucose concentration was increased and simultaneously GP and PKA activated. Apart from the slight dip in cAMP observed at time 10, the simulation accurately reproduces published experimental data [9, 19, 38].

Figure 5
figure 5

Predictions of Simplified cAMP Model A. (A) Species concentrations before and after a pulse of glucose. (B) Cyclic AMP levels of pde1Δ and pde2Δ mutants: blue trace – wild type; red trace – pde2Δ; green trace – pde1Δ. Glucose is increased to 5 after 5 seconds in both simulations.

To test if the model would also accurately reproduce phenotypic cAMP profiles of pde1Δ and pde2Δ mutants, the cAMP ODE (equation 10 defined above) was modified to remove the Pde1 and Pde2 reactions. The resultant "mutant" models were simulated as before, and as shown in Figure 5 (panel B), the simulations accurately reproduce the experimental data [19] (again with the exception of the slight dip in cAMP profile seen in the wild type and pde2Δ model mutants). We therefore conclude that this greatly simplified conceptual model is capable of reproducing the essential dynamics of changes in cAMP levels observed in response to glucose addition in wild-type as well as in the cAMP phosphodiesterase deletion mutants.

We then used this model to test the roles of the Krh proteins, which according to Harashima and Heitman [22] act by stabilizing the Ira proteins, whereas according to Peeters et al. [7, 45] they directly inhibit PKA. Initially we incorporated the Krh proteins into the Simplified cAMP Model A based on their function proposed by Peeters et al. The model extended in this way is called Simplified cAMP Model B. The ODE for PKA has been modified accordingly to include the Krh proteins (equation 13). The rate of PKA activation is decreased by Krh, and the rate of PKA deactivation is increased by Krh:

(13)

The formation of the G protein complex was modeled with the following mass action equation:

(14)

We tested Simplified cAMP Model B to see if it could reproduce the results from studies on adenylate cyclase mutants by Peeters et al. [7]. For this purpose, adenylate cyclase was removed from the model. The adenylate cyclase deletion model (cyr1Δ) was simulated with cAMP concentration set to 1. The GPA2Q 300L(Gpa2 constitutively active) mutant was modeled by setting the concentration of GPa to 1 and the parameter V maxGPdeact to 0. The pde2Δ mutant model was simulated as described earlier. The krh mutant was simulated fixing GPa levels to 1.

As shown in Figure 6 the cAMP and PKAa levels of the mutant model simulations are in agreement with experimentally observed phenotypes. The cyr1Δ model mutant has near-zero steady state levels of cAMP and PKA, which corresponds well with the fact that a cyr1Δ mutant is nonviable. Deleting Pde2p in the model elevates cAMP and PKAa levels, a result which agrees with the observation that a cyr1Δpde2Δ mutant is viable if supplemented with external cAMP. Deletion of Krh in the model produces a further increase in PKAa, which is in agreement with the observation that these mutants require less exogenous cAMP for viability [45]. Simulation of the model gives results that correspond well to the observations of Peeters et al., [7] when Krh is modeled as a direct inhibitor of PKA.

Figure 6
figure 6

Cyclic AMP and PKAa levels in Simplified cAMP Model B mutants when cAMP levels are set to 1 and PKAa set to 0. Model mutant genotypes are: cyr1Δ (green), cyr1Δpde2Δ (blue), cyr1Δpde2Δkrh1/2Δ (red), cyr1Δpde2Δ GPA2Q 300L(cyan).

We attempted to make a model of Krh activity as proposed by Harashima and Heitman [22]. In the Simplified cAMP Model B, Krh is quickly reassociated with the G proteins, allowing the system to exert negative feedback. However, any feedback in the mechanism proposed by Harashima and Heitman [22] is impossible because the Ira proteins are degraded, and re-synthesis of these proteins could not be fast enough to allow the Ira proteins to inhibit the Ras proteins. Therefore in all further developments of the complete cAMP pathway models Krh was retained as a direct inhibitor of PKA.

Modelling the complete cAMP pathway's response to glucose

The simplified conceptual model allowed us to capture the essential dynamics of the cAMP pathway and also to successfully incorporate the role of the Krh proteins. In order to fully understand the dynamics of this complex pathway we created a deterministic model which includes all components of the cAMP pathway and their interrelationships as currently reported. This Complete cAMP Model (Table 1 and Figure 7) consists of several distinct "modules" which are described below.

Figure 7
figure 7

Schematic representation of the Complete cAMP model, using the SBGN notation.

Glucose import and metabolism

The import of glucose was modelled using the following equation, as in [46]:

(15)

where v tr is the rate of transport (in mM per second), s is the extracellular glucose concentration, p is the intracellular glucose concentration, K M is the Michaelis constant (in mM) and K i is the interaction constant.

The metabolism of glucose via glycolysis was summarized with mass action kinetics, so that the intracellular glucose concentration did not exceed 1.5 mM during simulation, as described [46].

Gpa2 and Krh

As described earlier, Gpa2 is activated by Gpr1, and Gpr1 is activated by extracellular glucose. The activation of Gpr1 is modelled with mass action kinetics, whereby Gpr1 forms a complex with extracellular glucose. The activation of Gpa2 is based on mass action kinetics, with activated Gpr1 as an essential activator. Deactivation of Gpa2 is modelled using a basal rate of deactivation (representing the intrinsic GTPase activity of Gpa2), which can be enhanced by Rgs2. The binding of Gpa2 to Krh to form a complex is represented with simple mass action kinetics.

Ras2

Ras2 is very challenging to model because a large number of molecular species are involved in its regulation. It is directly activated by Cdc25, but it is activated indirectly by glucose. We chose to model the activation of Ras2 using general hyperbolic modifier kinetics. In this reaction, glucose acts as a modifier which increases the rate of the reaction, but the reaction is dependent on Cdc25. The deactivation of Ras2 was modelled using modified mass action kinetics with Ira as an activator. This captured the intrinsic GTPase activity of Ras2.

Adenylate Cyclase

Adenylate cyclase is represented as a Michaelis-Menten enzyme, with the following modifications. Activated Gpa2 and activated Ras increase the k cat of adenylate cyclase, increasing the maximum activity of the enzyme. In order to simplify the model, the substrate for adenylate cyclase (ATP) is not included, as the intracellular concentration of ATP is always greatly in excess of the cAMP concentration.

PKA

PKA is modelled using the mass action kinetics module with the addition of the actions of the Krh proteins described earlier. The forward reaction (PKA dissociation) is inhibited by Krh, and the backward reaction (PKA association) is activated by Krh.

The Phosphodiesterases

Pde2 is represented as a Michaelis Menten enzyme with a K m value of 0.002 mM, determined by parameter estimation (Table 4). Pde1 has been shown to be activated by phosphorylation, so the phosphorylated form has a lower K m and higher k cat than the dephosphorylated form. For this reason, Pde1 is represented by two species – the phosphorylated and the dephosphorylated form of Pde1p, respectively.

Table 4 Parameters of the complete cAMP pathway Model

The model was written in SBML format [36] and is included as Additional file 1. The pathway diagram was constructed using CellDesigner [47], incorporating the Systems Biology Graphical Notation (SBGN) scheme http://www.sbgn.org/. The representation of the model is shown in Figure 7, and its reactions and rate laws are shown in Table 5. The parameters of the Complete cAMP Model are given in Table 4, including both estimated and experimentally derived parameters. The cAMP data used for the parameter estimation were taken from [38], where 5 mM glucose was added to glucose-starved cell suspension after 60 seconds, followed by the addition of 100 mM glucose after 240 seconds. The cAMP profile (Figure 8) computed by simulation of our Complete cAMP Model after parameter estimation is in good agreement with previous observations.

Table 5 Reactions of the Complete cAMP pathway Model
Figure 8
figure 8

Result of parameter estimation of the Complete cAMP Model. The blue dotted trace represents cAMP data from [35], and the solid blue trace represents cAMP levels computed by the model.

The Complete cAMP Model illustrates several important features of the pathway. The balance of flux between cAMP synthesis and hydrolysis (Figure 9, panel A) demonstrates that Pde1p is more important than is Pde2p for controlling the cAMP levels following glucose pulses, as the effect of Pde1p on the rate of cAMP hydrolysis is much greater than that of Pde2p. Furthermore, the level of active Gpa2 is proportional to the level of extracellular glucose, and the level of Krh drops as it forms a complex with activated Gpa2 (Figure 9, panel B). Importantly, the proportion of active PKA (Figure 9, panel C) is not directly proportional to the cAMP level (Figure 8), allowing PKA to exert negative feedback on the cAMP level, even when the cAMP level drops. PKA exerts this feedback by activating Pde1p (Figure 9, panel A) and deactivating Ras2 via phosphorylation of Cdc25 (Figure 9, panel D). The latter mechanism is a feature of the Complete cAMP Model but not of the Simplified cAMP Model A and explains why cAMP level can come down after a glucose pulse in the Complete cAMP Model but not in the Simplified cAMP Model A.

Figure 9
figure 9

Predictions of the Complete cAMP Model. (A) cAMP synthesis and hydrolysis rates. Glucose is increased to 5 mM at time 60, and increased to 100 mM at time 240. Blue trace: rate of cAMP synthesis. Red trace: rate of cAMP hydrolysis by Pde1. Green trace: rate of cAMP hydrolysis by Pde2. (B) Levels of species in the Gpa2 module. Blue trace: inactive Gpa2. Green trace: active Gpa2. Red trace: Krh. Cyan trace: complex of activated Gpa2 and Krh. (C) Levels of active (blue trace) and inactive (green trace) PKA. (D) Levels of Ras2a and Cdc25 in response to 5 mM glucose at time 60sec, and 100 mM at 240 sec. Top part: Ras2a (blue trace); bottom part: phosphorylated (green trace) and unphosphorylated (blue trace) Cdc25.

Discussion

We have successfully created a series of deterministic mathematical models to investigate the cAMP pathway in S. cerevisiae. These range from simplified, conceptual models of the pathway, to an extensive model that fits experimental data. We were able to build a simplified model of the PKA module, containing only two variables and two parameters, without compromising the behaviour of the system. The simplification of the PKA module demonstrates the power of deterministic models. The components of this pathway are present in high abundance (proteins in thousands, nucleotides in millions per cell), making a deterministic model better suited than a stochastic one (we note also that we are not seeking to model potentially hundreds of kinds of protein molecule with different post-translational modifications).

In our PKA Model, the activation of PKA is worthy of particular attention. In previously published models, PKA activity was directly proportional to the cAMP level [14]. However, it has been proposed that PKA autophosphorylation provides a feed-forward mechanism for PKA activation [48], as Tpk1p is phosphorylated following a glucose pulse [21]. Alternatively, it is proposed that Krh inhibits PKA, and this inhibition is removed when Krh is recruited to activated Gpa2 [45]. Our Simplified cAMP Model B shows that the latter scenario is more likely, as this model corresponds well with observable phenotypes.

Our Simplified cAMP Model shows that the basic dynamics of the pathway in response to glucose can be explained with a relatively straightforward feedback mechanism. The activation of PKA by cAMP, followed by the activation of Pde1 and the inhibition of adenylate cyclase is sufficient to produce a characteristic "spike" of cAMP, followed by the emergence of a new steady state level of cAMP and PKA. This model has been tested by creating phosphodiesterase deletion mutant models (Figure 5, panel B). Deleting Pde2 in the model results in a higher steady state level of cAMP, but it does not significantly affect the cAMP spike. This phenotype is indeed found in yeast pde2Δ mutants [3]. However, removing Pde1 from the model results in a cAMP spike with increased peak height and duration, which is comparable to that experimentally determined in pde1Δ mutant [19].

In the Simplified cAMP model A, a slight dip in the level of cAMP can be seen before the cAMP level reaches a steady state after a pulse of glucose. Although this slight oscillation is not widely noted in the literature, it is possible to observe it in some experiments [38]. The presence of the slight oscillation in the model is dependent on the parameters of the model and the glucose concentration. It remains to be seen whether this oscillation is truly present in all or any circumstances.

The Simplified cAMP Model B (which incorporates the Krh proteins) demonstrates the significance of the negative feedback. Furthermore, it shows that this feedback is possible if the Krh proteins were acting as direct inhibitors of PKA as proposed by Peeters et al. [7, 45] rather than stabilising the Ira proteins as proposed by Harashima and Heitman [22]. At the same time, it predicts that cAMP levels should decrease more rapidly in the cyr1Δpde2Δkrh1Δkrh2Δ mutant than in the cyr1Δpde2Δ mutant. It will be interesting to see if these mutants behave in the way predicted by our models.

Although the Simplified cAMP Model could account for the majority of the behaviour of the cAMP pathway, there were exceptions. Most notably, in simulations of the pde1Δ mutant model, the steady state level of cAMP became significantly higher after a glucose pulse than it was before (Figure 5, panel B). This is not seen experimentally, where there is little difference between the post-glucose cAMP levels seen in a wild type and pde1Δ mutant [19]. This feature of the Simplified cAMP Model prompted us to develop the Complete Model. Our Complete cAMP Model represents the first effort to consolidate all the known elements of the cAMP pathway into one deterministic mathematical model. In addition to this, we have fitted the parameters of our model to experimental data. The fact that the complete cAMP pathway model can reproduce cAMP levels found in the literature indicates that the model is a reliable in silico approximation of the in vivo system. Furthermore, our model has other advantages. Firstly, as a deterministic model, it is computationally inexpensive to simulate and easy to analyze. Secondly, it represents a physiologically realistic steady state before glucose is introduced, in that the cAMP level is not zero. This contrasts with the model found in [14], in which the cAMP level is set to zero before glucose addition, which is biologically impossible, as cAMP is required for cell viability. After glucose addition, the model correctly represents the dynamical changes in cAMP level, until the cAMP level reaches a new steady state.

The models of the cAMP pathway described in this study make a number of predictions that could be tested experimentally. As a matter of further investigations in our lab, different species would be characterized following a pulse of glucose in terms of phosphorylation (Pde1p, Cdc25p), GTP loading (for Gpa2p), changes in cAMP levels (in cyr1Δpde2Δkrh1Δkrh2Δ in comparison to cyr1Δpde2Δ). Our Complete cAMP Model will no doubt be improved and tested further in the future. As more parameters are derived through experimentation, they can be included into the model to replace currently estimated parameters. We provide this model in SBML (Additional file 1), so that it can be easily expanded as scientific knowledge increases. For example, details on the mechanism of glucose activation of Ras2 could be incorporated when this mechanism is elucidated.

This model could be integrated with models of other pathways, a good example being that of the cell cycle, given the fact that cell cycle progression is controlled partly by the cAMP pathway [49]. It could also be integrated with a metabolic model such as the community consensus version recently published [50] via known PKA targets. Furthermore it could be adapted to other organisms such as the human fungal pathogen Candida albicans, as it is well documented that the cAMP pathway plays a key role in regulating virulence [51].

Conclusion

We report a deterministic mathematical model of the cAMP-mediated signal transduction pathway in S. cerevisiae. The model is easier to compute and simulate as it has a reduced number of variables and parameters in comparison to previously reported stochastic model of this pathway. Furthermore, our model contains components such as the regulatory Krh proteins that have not been included before. It is able to simulate accurately experimentally derived patterns of cAMP changes observed in different pathway mutants in response to glucose addition. We suggest that it is suitable for integration with other models such as that of the cell cycle or metabolism and that it could be adapted to medically important yeast species such as the human fungal opportunistic pathogen C. albicans.

References

  1. Portela P, Howell S, Moreno S, Rossi S: In vivo and in vitro phosphorylation of two isoforms of yeast pyruvate kinase by protein kinase A. J Biol Chem. 2002, 277: 30477-30487.

    Article  CAS  PubMed  Google Scholar 

  2. Longo VD: he Ras and Sch9 pathways regulate stress resistance and longevity. Exp Gerontol. 2003, 38: 807-811.

    Article  CAS  PubMed  Google Scholar 

  3. Park J, Grant CM, Dawes IW: The high-affinity cAMP phosphodiesterase of Saccharomyces cerevisiae is the major determinant of cAMP levels in stationary phase: involvement of different branches of the Ras-cyclic AMP pathway in stress responses. Biochem Biophys Res Commun. 2005, 327: 311-319.

    Article  CAS  PubMed  Google Scholar 

  4. Schneper L, Krauss A, Miyamoto R, Fang S, Broach JR: The Ras/protein kinase A pathway acts in parallel with the Mob2/Cbk1 pathway to effect cell cycle progression and proper bud site selection. Eukaryot Cell. 2004, 3: 108-120.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  5. Pedruzzi I, Dubouloz F, Cameroni E, Wanke V, Roosen J, Winderickx J, De Virgilio C: TOR and PKA signaling pathways converge on the protein kinase Rim15 to control entry into G0. Mol Cell. 2003, 12: 1607-1613.

    Article  CAS  PubMed  Google Scholar 

  6. Thevelein JM, Geladé R, Holsbeeks I, Lagatie O, Popova Y, Rolland F, Stolz F, Velde Van de S, Van Dijck P, Vandormael P, Van Nuland A, Van Roey K, Van Zeebroeck G, Yan B: Nutrient sensing systems for rapid activation of the protein kinase A pathway in yeast. Biochem Soc Trans. 2005, 33 (Pt 1): 253-256.

    Article  CAS  PubMed  Google Scholar 

  7. Peeters T, Louwet W, Geladé R, Nauwelaers D, Thevelein JM, Versele M: Kelch-repeat proteins interacting with the Galpha protein Gpa2 bypass adenylate cyclase for direct regulation of protein kinase A in yeast. Proc Natl Acad Sci USA. 2006, 103: 13034-13039.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  8. Nikawa J, Cameron S, Toda T, Ferguson KM, Wigler M: Rigorous feedback control of cAMP levels in Saccharomyces cerevisiae. Genes Dev. 1987, 1: 931-937.

    Article  CAS  PubMed  Google Scholar 

  9. Kraakman L, Lemaire K, Ma P, Teunissen AW, Donaton MC, Van Dijck P, Winderickx J, de Winde JH, Thevelein JM: A Saccharomyces cerevisiae G-protein coupled receptor, Gpr1, is specifically required for glucose activation of the cAMP pathway during the transition to growth on glucose. Mol Microbiol. 1999, 32: 1002-1012.

    Article  CAS  PubMed  Google Scholar 

  10. Versele M, de Winde JH, Thevelein JM: A novel regulator of G protein signalling in yeast, Rgs2, downregulates glucose-activation of the cAMP pathway through direct inhibition of Gpa2. EMBO J. 1999, 18: 5577-5591.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  11. Lai CC, Boguski M, Broek D, Powers S: Influence of guanine nucleotides on complex formation between Ras and CDC25 proteins. Mol Cell Biol. 1993, 13: 1345-1352.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  12. Crechet JB, Poullet P, Mistou MY, Parmeggiani A, Camonis J, Boy-Marcotte E, Damak F, Jacquet M: Enhancement of the GDP-GTP exchange of RAS proteins by the carboxyl-terminal domain of SDC25. Science. 1990, 248: 866-868.

    Article  CAS  PubMed  Google Scholar 

  13. Colombo S, Ronchetti D, Thevelein JM, Winderickx J, Martegani E: Activation state of the Ras2 protein and glucose-induced signaling in Saccharomyces cerevisiae. J Biol Chem. 2004, 279: 46715-46722.

    Article  CAS  PubMed  Google Scholar 

  14. Cazzaniga P, Pescini D, Besozzi D, Mauri G, Colombo S, Martegani E: Modeling and stochastic simulation of the Ras/cAMP/PKA pathway in the yeast Saccharomyces cerevisiae evidences a key regulatory function for intracellular guanine nucleotides pools. J Biotechnol. 2007, 133: 377-385.

    Article  PubMed  Google Scholar 

  15. Sass P, Field J, Nikawa J, Toda T, Wigler M: Cloning and characterization of the high-affinity cAMP phosphodiesterase of Saccharomyces cerevisiae. Proc Natl Acad Sci USA. 1986, 83: 9303-9307.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  16. Ishikawa T, Matsumoto K, Uno I: Yeast mutants altered in the cAMP cascade system. Methods Enzymol. 1988, 159: 27-42.

    Article  CAS  PubMed  Google Scholar 

  17. Uno I, Matsumoto K, Ishikawa T: Characterization of a cyclic nucleotide phosphodiesterase-deficient mutant in yeast. J Biol Chem. 1983, 258: 3539-3542.

    CAS  PubMed  Google Scholar 

  18. Ptacek J, Devgan G, Michaud G, Zhu H, Zhu X, Fasolo J, Guo H, Jona G, Breitkreutz A, Sopko R, McCartney RR, Schmidt MC, Rachidi N, Lee S, Mah AS, Meng L, Stark MJR, Stern DF, De Virgilio C, Tyers M, Andrews B, Gerstein M, Schweitzer B, Predki PF, Snyder M: Global analysis of protein phosphorylation in yeast. Nature. 2005, 438: 679-684.

    Article  CAS  PubMed  Google Scholar 

  19. Ma P, Wera S, Van Dijck P, Thevelein JM: The PDE1-encoded low-affinity phosphodiesterase in the yeast Saccharomyces cerevisiae has a specific function in controlling agonist-induced cAMP signaling. Mol Biol Cell. 1999, 10: 91-104.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  20. Gross E, Goldberg D, Levitzki A: Phosphorylation of the S. cerevisiae Cdc25 in response to glucose results in its dissociation from Ras. Nature. 1992, 360: 762-765.

    Article  CAS  PubMed  Google Scholar 

  21. Portela P, Moreno S: Glucose-dependent activation of protein kinase A activity in Saccharomyces cerevisiae and phosphorylation of its TPK1 catalytic subunit. Cell Signal. 2006, 18: 1072-1086.

    Article  CAS  PubMed  Google Scholar 

  22. Harashima T, Anderson S, Yates JR, Heitman J: The kelch proteins Gpb1 and Gpb2 inhibit Ras activity via association with the yeast RasGAP neurofibromin homologs Ira1 and Ira2. Mol Cell. 2006, 22: 819-830.

    Article  CAS  PubMed  Google Scholar 

  23. Lu A, Hirsch JP: Cyclic AMP-independent regulation of protein kinase A substrate phosphorylation by Kelch repeat proteins. Eukaryot Cell. 2005, 4: 1794-1800.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Kell DB, Oliver SG: Here is the evidence, now what is the hypothesis? The complementary roles of inductive and hypothesis-driven science in the post-genomic era. Bioessays. 2004, 26: 99-105.

    Article  PubMed  Google Scholar 

  25. Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI – a COmplex PAthway SImulator. Bioinformatics. 2006, 22: 3067-3074.

    Article  CAS  PubMed  Google Scholar 

  26. Kell DB, Knowles JD: The role of modeling in systemsbiology. System modeling in cellular biology: from concepts to nuts and bolts. Edited by: Szallazi Z, Stelling J, Periwal V. 2006, 3-18. Cambridge, MIT Press

    Chapter  Google Scholar 

  27. Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10: 122-133.

    Article  CAS  PubMed  Google Scholar 

  28. Teusink B, Passarge J, Reijenga CA, Esgalhado E, Weijden van der CC, Schepper M, Walsh MC, Bakker BM, van Dam K, Westerhoff HV, Snoep JL: Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry. Eur J Biochem. 2000, 267: 5313-5329.

    Article  CAS  PubMed  Google Scholar 

  29. Kofahl B, Klipp E: Modelling the dynamics of the yeast pheromone pathway. Yeast. 2004, 21: 831-850.

    Article  CAS  PubMed  Google Scholar 

  30. Yu RC, Pesce CG, Colman-Lerner A, Lok L, Pincus D, Serra E, Holl M, Benjamin K, Gordon A, Brent R: Negative feedback that improves information transmission in yeast signalling. Nature. 2008, 456: 755-761.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  31. Yu RC, Resnekov O, Abola AP, Andrews SS, Benjamin KR, Bruck J, Burbulis IE, Colman-Lerner A, Endy D, Gordon A, Holl M, Lok L, Pesce CG, Serra E, Smith RD, Thomson TM, Tsong AE, Brent R: The Alpha Project: a model system for systems biology research. IET Syst Biol. 2008, 2: 222-233.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  32. Csikasz-Nagy A, Battogtokh D, Chen KC, Novak B, Tyson JJ: Analysis of a generic model of eukaryotic cell-cycle regulation. Biophys J. 2006, 90: 4361-4379.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  33. Longo D, Hasty J: Dynamics of single-cell gene expression. Mol Syst Biol. 2006, 2: 64-

    Article  PubMed Central  PubMed  Google Scholar 

  34. Salwinski L, Eisenberg D: In silico simulation of biological network dynamics. Nat Biotechnol. 2004, 22: 1017-1019.

    Article  CAS  PubMed  Google Scholar 

  35. Mendes P: GEPASI: a software package for modelling the dynamics, steady states and control of biochemical and other systems. Comput Appl Biosci. 1993, 9: 563-571.

    CAS  PubMed  Google Scholar 

  36. Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, Cornish-Bowden A, Cuellar AA, Dronov S, Gilles ED, Ginkel M, Gor V, Goryanin II, Hedley WJ, Hodgman TC, Hofmeyr JH, Hunter PJ, Juty NS, Kasberger JL, Kremling A, Kummer U, Le Novere N, Loew LM, Lucio D, Mendes P, Minch E, Mjolsness ED, Nakayama Y, Nelson MR, Nielsen PF, Sakurada T, Schaff JC, Shapiro BE, Shimizu TS, Spence HD, Stelling J, Takahashi K, Tomita M, Wagner J, Wang J: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19: 524-531.

    Article  CAS  PubMed  Google Scholar 

  37. Schmidt H, Jirstrand M: Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology. Bioinformatics. 2006, 22: 514-515.

    Article  CAS  PubMed  Google Scholar 

  38. Rolland F, De Winde JH, Lemaire K, Boles E, Thevelein JM, Winderickx J: Glucose-induced cAMP signalling in yeast requires both a G-protein coupled receptor system for extracellular glucose detection and a separable hexose kinase-dependent sensing process. Mol Microbiol. 2000, 38: 348-358.

    Article  CAS  PubMed  Google Scholar 

  39. Davey HM, Kell DB: Flow cytometry and cell sorting of heterogeneous microbial populations: the importance of single-cell analyses. Microbiol Rev. 1996, 60: 641-696.

    PubMed Central  CAS  PubMed  Google Scholar 

  40. Cerny V: A thermodynamical approach to the travelling salesman problem: An efficient simulation algorithm. Journal of Optimization Theory and Applications. 1985, 45: 41-51.

    Article  Google Scholar 

  41. Kirkpatrick S, Gelatt CD, Vecchi MP: Optimization by Simulated Annealing. Science. 1983, 220 (4598): 671-680.

    Article  CAS  PubMed  Google Scholar 

  42. Danø S, Madsen MF, Schmidt H, Cedersund G: Reduction of a biochemical model with preservation of its basic dynamic properties. FEBS J. 2006, 273: 4862-4877.

    Article  PubMed  Google Scholar 

  43. Toda T, Cameron S, Sass P, Zoller M, Wigler M: Three different genes in S. cerevisiae encode the catalytic subunits of the cAMP-dependent protein kinase. Cell. 1987, 50: 277-287.

    Article  CAS  PubMed  Google Scholar 

  44. Matsumoto K, Uno I, Oshima Y, Ishikawa T: Isolation and characterization of yeast mutants deficient in adenylate cyclase and cAMP-dependent protein kinase. Proc Natl Acad Sci USA. 1982, 79: 2355-2359.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  45. Peeters T, Versele M, Thevelein JM: Directly from Galpha to protein kinase A: the kelch repeat protein bypass of adenylate cyclase. Trends Biochem Sci. 2007, 32: 547-554.

    Article  CAS  PubMed  Google Scholar 

  46. Teusink B, Diderich JA, Westerhoff HV, van Dam K, Walsh MC: Intracellular glucose concentration in derepressed yeast cells consuming glucose is high enough to reduce the glucose transport rate by 50%. J Bacteriol. 1998, 180: 556-562.

    PubMed Central  CAS  PubMed  Google Scholar 

  47. Funahashi A, Jouraku A, Matsuoka Y, Kitano H: Integration of CellDesigner and SABIO-RK. In Silico Biol. 2007, 7 (2 Suppl ): S81-S90.

    PubMed  Google Scholar 

  48. Vaseghi S, Macherhammer F, Zibek S, Reuss M: Signal transduction dynamics of the protein kinase-A/phosphofructokinase-2 system in Saccharomyces cerevisiae. Metab Eng. 2001, 3: 163-172.

    Article  CAS  PubMed  Google Scholar 

  49. Tyson JJ, Novak B: Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions. J Theor Biol. 2001, 210: 249-263.

    Article  CAS  PubMed  Google Scholar 

  50. Herrgård MJ, Swainston N, Dobson P, Dunn WB, Arga KY, Arvas M, Blüthgen N, Borger S, Costenoble R, Heinemann M, Hucka M, Le Novère N, Li P, Liebermeister W, Mo ML, Oliveira AP, Petranovic D, Pettifer S, Simeonidis E, Smallbone K, Spasiã I, Weichart D, Brent R, Broomhead DS, Westerhoff HV, Kirdar B, Penttilä M, Klipp E, Palsson BØ, Sauer U, Oliver SG, Mendes P, Nielsen J, Kell DB: A consensus yeast metabolic network reconstruction obtained from a community approach to systems biology. Nat Biotechnol. 2008, 26: 1155-1160.

    Article  PubMed Central  PubMed  Google Scholar 

  51. Wilson D, Tutulan-Cunita A, Jung WH, Hauser N, Henrnandez R, Williamson T, Piekarska K, Rupp S, Young T, Stateva L: Deletion of the high-affinity cAMP phosphodiesterase encoded by PDE2 affects stress responses and virulence in Candida albicans. Mol Microbiol. 2007, 65: 841-856.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

The authors acknowledge the helpful recommendations of the anonymous reviewers.

TW is a grateful recipient of a BBSRC-funded PhD studentship. JMS, LS and TW thank the BBSRC and EPSRC for support of the Manchester Centre for Integrative Systems Biology http://www.mcisb.org/.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Lubomira Stateva.

Additional information

Authors' contributions

TW carried out the modeling work and drafted the manuscript. J-MS and DK participated in different stages of the modeling work and in writing of the manuscript. LS conceived of the study and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Williamson, T., Schwartz, JM., Kell, D.B. et al. Deterministic mathematical models of the cAMP pathway in Saccharomyces cerevisiae. BMC Syst Biol 3, 70 (2009). https://0-doi-org.brum.beds.ac.uk/10.1186/1752-0509-3-70

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/1752-0509-3-70

Keywords