Parameter-free model discrimination criterion based on steady-state coplanarity

Heather A. Harrington$^{*}$, Kenneth L. Ho$^{\dagger}$, Thomas Thorne$^{*}$, Michael P. H. Stumpf$^{*}$

$^{*}$Theoretical Systems Biology, Division of Molecular Biosciences, Imperial College London, Wolfson Building, London, SW7 2AZ, UK

$^{\dagger}$Courant Institute of Mathematical Sciences and Program in Computational Biology, New York University, 251 Mercer Street, New York, NY 10012

About this worksheet

This worksheet contains all computations referenced in:

Harrington HA, Ho KL, Thorne T, Stumpf MPH (2012) Parameter-free model discrimination criterion based on steady-state coplanarity. Proc Natl Acad Sci USA 109(39): 15746–15751. doi:10.1073/pnas.1117073109.

which were performed using Sage version:

Please read the associated paper for details.

Sage was chosen for its strength in unifying many different software packages under a central frontend, providing capabilities for symbolic and numerical computation, both of which are used extensively here. Worksheet cells typically make use of previously computed results; cells are therefore meant to be evaluated in order of their appearance.

As a preliminary, we import all required packages and define some settings for plotting.

Multisite phosphorylation

We first consider the two-site phosphorylation model from [1]:

$$\begin{align} K + S_{u} &\mathop{\longleftrightarrow}^{a_{u}}_{b_{u}} K S_{u} \mathop{\longrightarrow}^{c_{uv}} K + S_{v},\\ F + S_{v} &\mathop{\longleftrightarrow}^{\alpha_{v}}_{\beta_{v}} F S_{v} \mathop{\longrightarrow}^{\gamma_{vu}} F + S_{u}, \end{align}$$

where $u, v \in \{ 0, 1 \}^{2}$ are bit strings of length two, encoding the occupancies of each site ($0$ or $1$ for the absence or presence, respectively, of a phosphate), with $u$ having less bits than $v$; $S_{u}$ is the phosphoform with phosphorylation state $u$; $K$ is a kinase, an enzyme that adds phosphates; and $F$ is a phosphatase, an enzyme that removes phosphates.

Each enzyme can be either processive (P), where more than one phosphate modification may be achieved in a single step, or distributive (D), where only one modification is allowed before the enzyme dissociates from the substrate ($c_{0011} = 0$ for $K$, $\gamma_{1100} = 0$ for $F$). This mechanistic diversity generates four competing models: PP, PD, DP, and DD; where the first letter designates the mechanism of the kinase, and the second, that of the phosphatase.

Below, we define the required variables (throughout, we follow the convention that lowercase letters denote the concentrations of the corresponding species indicated in uppercase) and parameters

$$\begin{align} \mathbf{x} &= \left( k s_{00}, k s_{01}, k s_{10}, f s_{01}, f s_{10}, f s_{11}, k, f, s_{00}, s_{01}, s_{10}, s_{11} \right),\\ \mathbf{a} &= \left( a_{00}, a_{01}, a_{10}, b_{00}, b_{01}, b_{10}, c_{0001}, c_{0010}, c_{0011}, c_{0111}, c_{1011}, \alpha_{01}, \alpha_{10}, \alpha_{11}, \beta_{01}, \beta_{10}, \beta_{11}, \gamma_{0100}, \gamma_{1000}, \gamma_{1100}, \gamma_{1101}, \gamma_{1110} \right), \end{align}$$

respectively, where we have ordered $\mathbf{x}$ in anticipation of eliminating all variables except the $s_{u}$ using Gröbner basis methods.

We construct the polynomial ring $\mathbb{Q} [\mathbf{a}]$ in the parameters $\mathbf{a}$ with coefficients from the rational numbers $\mathbb{Q}$. The field over which we will perform all Gröbner basis calculations is the fraction field $\mathbb{K}$ of $\mathbb{Q} [\mathbf{a}]$, i.e., the polynomials generating the corresponding ideal are elements of $\mathbb{K} [\mathbf{x}]$.

We then define the dynamics $\dot{\mathbf{x}} = \mathbf{f} (\mathbf{x})$ governing the PP system, assuming mass-action kinetics. The dynamics of the other models will be derived from this by setting $c_{0011} = 0$ or $\gamma_{1100} = 0$ as appropriate.

As a check, we make sure that the following conservations are satisfied:

To obtain the dynamics for the PD, DP, and DD models, we need to apply certain replacement rules to the dynamics of the PP model. Sage has an internal facility for assigning variables to specific values, but there is some issue when the base rings are complicated. Below, we define a small hack that allows this to be done.

We now compute the reduced Gröbner basis of each model with respect to the lexicographic ordering over the field $\mathbb{K}$.

The PP model is first. There are twelve polynomials in the Gröbner basis, with variables:

Next is the PD model, with variables:

Then comes the DP model, with variables:

Finally, we consider the DD model, whose basis polynomials have variables:

Each model evidently has twelve basis polynomial, none of which contain only the variables $s_{u}$. The last three come close, but all contain also the phosphatase concentration $f$. As in [1], however, we find that all three polynomials are of the form $p(f, s_{00}, s_{01}, s_{10}, s_{11}) = f \cdot q(s_{00}, s_{01}, s_{10}, s_{11})$:

hence implying that we can simply divide through by $f \neq 0$ ($f = 0$ only if there is no phosphatase in the system, which we assume not to be the case). For convenience, we also multiply through by the denominator so that all coefficients are elements of $\mathbb{Q} [\mathbf{a}]$ instead of $\mathbb{K}$ (this avoids "division by zero" issues in invariant minimization, though we have to check that we are not multiplying by zero and trivializing the invariant). This gives us three steady-state invariants for each model, in the variables $s_{00}$, $s_{01}$, $s_{10}$, and $s_{11}$ only.

We check that the resulting polynomials consist only of the variables $s_{00}$, $s_{01}$, $s_{10}$, and $s_{11}$ by listing the monomials in each invariant:

Now that we have our steady-state invariants, we need to generate data on which to test our model selection heuristics. We do this by integrating the dynamics for each model above until steady state, using the solver LSODA from ODEPACK. Below, we define the dynamics functions required for LSODA, using Cython for acceleration.

First, the PP dynamics:

Then the PD dynamics:

Next, the DP dynamics:

And, finally, the DD dynamics:

We generate random parameters for the models, sampled from a log-normal distribution with median $\mu^{*} = e^{\mu} = 1$ and multiplicative standard deviation $\sigma^{*} = e^{\sigma} = 2$, where $\mu$ and $\sigma$ are the mean and standard deviation, respectively, of the underlying normal distribution. Using these parameters, we generate $m = 100$ time course trajectories for each model by integration of the model ODEs over the time interval $0 \leq t \leq 100$; each trajectory is seeded with random initial conditions, also sampled from a log-normal distribution with $\mu^{*} = 1$ and $\sigma^{*} = 2$.

In anticipation of our experiments with noise, we also initialize our relative errors below.

We now integrate the four models.

To simulate measurement error, we corrupt the data with various levels $\epsilon$ of noise by drawing random multiplicative factors from log-normal distributions with $\mu^{*} = 1$ and $\sigma^{*} = 1 + \epsilon$, i.e., $\epsilon$ is the relative (multiplicative) error.

We save the data for reference.

To facilitate the analysis, we load the precomputed data. This is simply for the purpose of saving time; further simulations have indicated that the saved data are representative of the general case.

To compute the coplanarity error $\Delta$, we need to consider the Jacobian of the variable transformation $\varphi$ taking the original variables to the transformed variables in each steady-state invariant. This is implemented below along with various auxiliary functions.

We now compute the coplanarity error $\Delta$ for each model over each set of time course trajectories. The required computation is essentially a singular value decomposition over the data as organized by the steady-state invariants.

The computation is quite lengthy, so we separate the process into four parts. First, we analyze each invariant on data from the PP model:

Then we do the same for data from the PD model:

Next for the DP model:

And finally for the DD model:

As before, we save the coplanarity errors for reference.

For subsequent analysis, we load the precomputed data.

The coplanarity errors for each model/data pair are shown below at $\epsilon = 10^{-9}$, from which it is clear that the DP/DD models can be rejected on the basis of PP/PD data.

For comparison, the $\Delta$ cutoff at a significance level of $\alpha = 0.05$ is:

We also plot the evolution of the coplanarity errors for each model pair on PP data as a function of time. The intuition is that although coplanarity is a steady-state condition, some approximation of coplanarity should be evident provided that the data are reasonably close to steady state. The results suggest that the method can also be effective for time course data near steady state, provided that there is not too much noise.

The results clearly demonstrate that our rejection power decreases with noise. To estimate the critical noise level below which the DP/DD models can be rejected from PP data, we plot the DP/DD coplanarity errors at steady state as a function of $\epsilon$, which suggests a cutoff of $\epsilon \sim 10^{-4}$ for $\alpha = 0.05$.

Next, we try using invariant minimization as a discrimination tool between the phosphorylation models. We hence define the invariant error and its log-likelihood for each model below, consisting of the values of all invariants over all data, with each value scaled as appropriate:

First, for the PP model:

Then for the PD model:

Next for the DP model:

And, finally, for the DD model:

We now maximize the invariant error likelihood for each model over the model parameters using L-BFGS-B. This maximization proceeds in two phases. First, we obtain an approximate optimal parameter set by minimizing the normalized invariant error, using an initial guess of $1$ for each variable. This is then taken as an initial parameter estimate for the likelihood maximization. Lower and upper bounds of $0.01$ and $100$, respectively, were applied for all variables in each optimization.

As with the coplanarity error, we separate the computation into four parts. First, we analyze data from the PP model:

Then data from the PD model:

Next, data from the DP model:

And, finally, data from the DD model:

Again, we save the data for reference:

and, again, we load the precomputed data for ease of analysis.

The results of the optimization, which are essentially the Akaike information criteria (AIC) for each model, are shown below at $\epsilon = 10^{-9}$, from which we see that the correct model can be identified from the data. Recall that models with a lower AIC are preferred.

The following plot shows the invariant error AIC for each model/data pair as a function of noise, which suggests a critical noise level of $\epsilon \sim 10^{-5}$.

Cell death signaling

We now turn to our second example, which involves receptor-mediated cell death signaling, the so-called extrinsic pathway of apoptosis. Specifically, we focus on the formation of the death-inducing signaling complex (DISC), a multi-protein oligomer composed of FasL, a death ligand, and its cognate receptor Fas, for which we study two competing models.

The first [2], which we call the crosslinking model, is based on the successive binding of Fas ($R$) to FasL ($L$):

$$\begin{align} L + R &\mathop{\longleftrightarrow}^{3 k_{f}}_{k_{r}} C_{1},\\ C_{1} + R &\mathop{\longleftrightarrow}^{2 k_{f}}_{2 k_{r}} C_{2},\\ C_{2} + R &\mathop{\longleftrightarrow}^{k_{f}}_{3 k_{r}} C_{3}, \end{align}$$

where $C_{i}$ is the complex FasL:Fas$_{i}$.

The second [3], which we call the cluster model, posits three forms of Fas (inactive, $X$; active and unstable, $Y$; active and stable, $Z$) and specifies receptor cluster-stabilization events driven by FasL:

$$\begin{align} X &\mathop{\longleftrightarrow}^{k_{o}}_{k_{u}} Y,\\ Z &\mathop{\longrightarrow}^{k_{u}} Y,\\ jY + \left( i - j \right) Z &\mathop{\longrightarrow}^{k_{s}^{\left( i \right)}} \left( j - k \right) Y + \left( i - j + k \right) Z,\\ L + jY + \left( i - j \right) Z &\mathop{\longrightarrow}^{k_{l}^{\left( i \right)}} L + \left( j - k \right) Y + \left( i - j + k \right) Z, \end{align}$$

where the last two reactions represent entire families generated by taking $i = 2$ or $3$, with $j = 1, \dots, i$ and $k = 1, \dots, j$.

We compute the dynamics, Gröbner bases, and steady-state invariants for each model below. For steady-state invariants, we identify corresponding variables by considering the apoptotic signal $\zeta$ transduced by the DISC, defined as $\zeta = c_{1} + 2 c_{2} + 3 c_{3}$ for the crosslinking model, and $\zeta = z$ for the cluster model; the total FasL concentration ($\lambda = l + c_{1} + c_{2} + c_{3}$ and $\lambda = l$ for the crosslinking and cluster models, respectively); and the total Fas concentration ($\rho = r + c_{1} + 2 c_{2} + 3 c_{3}$ and $\rho = x + y + z$, respectively).

In all following code, $\lambda$ is represented as $\Lambda$ since the keyword lambda is reserved.

We begin with the crosslinking model. As for phosphorylation, we define the model variables and parameters, and the field over which Gröbner bases will be computed.

The dynamics of the model are given below:

We check for conservation of ligand and receptor:

Then we redefine the system dynamics in terms of $\lambda$, $\rho$, and $\zeta$, and compute a Gröbner basis using the ordering $(c_{2}, c_{3}, \lambda, \rho, \zeta)$.

Clearly, the fourth basis polynomial is suitable; we call that our steady-state invariant.

We now repeat the same process for the cluster model, first defining the model variables and parameters, and the base field.

The system dynamics are given as follows:

We check for conservation of Fas (FasL is conserved trivially):

We then compute a Gröbner basis for the ideal generated by the system dynamics using the variable ordering $(y, \lambda, \rho, \zeta)$.

As with the crosslinking model, it is immediate that only one steady-state invariant exists.

We list the monomials in each invariant below:

To generate data from the models, we define first the crosslinking model dynamics, implemented in Cython:

and then the cluster model dynamics:

For each model, we generate parameters from a log-normal distribution with $\mu^{*} = 1$ and $\sigma^{*} = 2$. Again, we generate $m = 100$ trajectories for each model, each with initial conditions also drawn from a log-normal distribution with $\mu^{*} = 1$ and $\sigma^{*} = 2$, and again we corrupt the resulting data with various levels $\epsilon$ of noise.

We integrate the model equations below. Some post-processing must be performed to obtain the values of the corresponding variables $\lambda$, $\rho$, and $\zeta$.

We then corrupt the data with noise.

We save the generated data for reference.

For ease of analysis, we load the precomputed data. Additional simulations show that these data, too, are representative of the general case.

The coplanarity error for each model/data time course pair is computed next:

We save the data here:

and, again, load it for subsequent analysis.

The results at $\epsilon = 10^{-9}$ show that the crosslinking model can be rejected on the basis of data from the cluster model.

We plot the time evolution of the coplanarity error $\Delta$ for each model on cluster model data below at various noise levels.

As with phosphorylation, we also study the behavior of $\Delta$ on cluster model data at steady state as a function of $\epsilon$. The critical noise level here is found to be $\epsilon \sim 10^{-2}$.

Finally, to test invariant minimization, we define the invariant error functions for each model. First, the crosslinking model:

and, then, the cluster model:

Invariant minimization is now performed, using the same protocol as for phosphorylation.

The data are saved for reference:

and loaded for subsequent analysis.

The AIC results below at $\epsilon = 10^{-9}$ show that the correct model can be identified easily from the data.

We also plot the invariant error AIC for each model/data pair below, which suggests that invariant minimization has discriminative power below $\epsilon \sim 10^{-2}$.

References

1. Manrai AK, Gunawardena J (2008) The geometry of multisite phosphorylation. Biophys J 95:5533–5543.
2. Lai R, Jackson TL (2004) A mathematical model of receptor-mediated apoptosis: dying to know why FasL is a trimer. Math Biosci Eng 1:325–338.
3. Ho KL, Harrington HA (2010) Bistability in apoptosis by receptor clustering. PLoS Comput Biol 6:e1000956.