Department of Mathematics and Centre for Integrative Systems Biology at Imperial College

About this worksheet

This worksheet takes advantage of the symbolic focus of Sage, additionally using PyLab (NumPy/SciPy + matplotlib) for numerical computation and plotting. Some worksheet cells make use of previously computed results; cells are therefore meant to be evaluated in order of their appearance.

Model formulation and analysis

Label each cluster by the tuple $(L, X, Y, Z)$, where $L$ is FasL, and $X$, $Y$, and $Z$ are three posited forms of Fas, denoting closed, open and unstable, and open and stable receptors, respectively. Within each cluster, define the reactions: $$\begin{align*} X &\mathop{\longleftrightarrow}^{k_{o}}_{k_{c}} Y\\ Z &\xrightarrow{k_{u}} Y,\\ jY + \left( i - j \right) Z &\xrightarrow{k_{s}^{\left( i \right)}} \left( j - k \right) Y + \left( i - j + k \right) Z, \quad \begin{cases} i = 2, \dots, m,\\ j = 1, \dots, i,\\ k = 1, \dots, j \end{cases}\\ L + jY + \left( i - j \right) Z &\xrightarrow{k_{l}^{\left( i \right)}} L + \left( j - k \right) Y + \left( i - j + k \right) Z, \quad \begin{cases} i = 2, \dots, n,\\ j = 1, \dots, i,\\ k = 1, \dots, j. \end{cases} \end{align*}$$

With the nondimensionalizations $$\xi = \frac{x}{s}, \quad \eta = \frac{y}{s}, \quad \zeta = \frac{z}{s}, \quad \lambda = \frac{l}{s}, \quad \tau = k_{c} t,$$ following the convention that lowercase letters denote the concentrations of their uppercase counterparts, and where $s$ is a characteristic concentration and $t$ is time, and $$\kappa_{o} = \frac{k_{o}}{k_{c}}, \quad \kappa_{u} = \frac{k_{u}}{k_{c}}, \quad \kappa_{s}^{\left( i \right)} = \frac{k_{s}^{\left( i \right)} s^{i - 1}}{k_{c}}, \quad \kappa_{l}^{\left( i \right)} = \frac{k_{l}^{\left( i \right)} s^{i}}{k_{c}},$$ the mass-action dynamics are $$\begin{align*} \frac{d \xi}{d \tau} &= -\nu_{o},\\ \frac{d \eta}{d \tau} &= \nu_{o} + \nu_{u} - \nu_{s} - \nu_{l},\\ \frac{d \zeta}{d \tau} &= -\nu_{u} + \nu_{s} + \nu_{l} \end{align*}$$ where $$\begin{align*} \nu_{o} &= \kappa_{o} \xi - \eta,\\ \nu_{u} &= \kappa_{u} \zeta,\\ \nu_{s} &= \sum_{i = 2}^{m} \kappa_{s}^{\left( i \right)} \sum_{j = 1}^{i} \eta^{j} \zeta^{i - j} \sum_{k = 1}^{j} k,\\ \nu_{l} &= \lambda \sum_{i = 2}^{n} \kappa_{l}^{\left( i \right)} \sum_{j = 1}^{i} \eta^{j} \zeta^{i - j} \sum_{k = 1}^{j} k. \end{align*}$$ Therefore, at steady state, $$\xi_{\infty} = \frac{\sigma - \zeta_{\infty}}{1 + \kappa_{o}}, \quad \eta_{\infty} = \kappa_{o} \xi_{\infty},$$ where $\sigma = \xi + \eta + \zeta$, and $\zeta_{\infty}$ is given by solving $d \zeta / d \tau = 0$ with $(\xi, \eta) \mapsto (\xi_{\infty}, \eta_{\infty})$, a polynomial in $\zeta_{\infty}$ of degree $\max \{ m, n \}$.

We then set baseline parameter values, and use Sage to perform various symbolic computations that will be useful later.

Our first figure shows the behavior of the model with $\lambda$ at various values of $\sigma$. The steady-state curves show bistability and hysteresis, with irreversible bistability at high $\sigma$. The stability of each steady state is computed by evaluating $(\partial / \partial \zeta) d \zeta / d \tau$ at $\zeta = \zeta_{\infty}$ (stable, solid lines; unstable, dashed lines).

The steady-state surface given as a parameterization over $(\lambda, \sigma)$ is shown below (blue, stable; red, unstable).

To further study the qualitative structure of the model, we look for the boundaries of the bistable regime in $(\lambda, \sigma)$-space. We approach this by considering the boundaries in $\sigma$ at a given value of $\lambda$, which we compute using a binary search algorithm. The two monostable regimes are characterized by the value of $\partial^{2} \lambda / \partial \zeta_{\infty}^{2}$, which is negative on the lower (life) branch and positive on the upper (death) branch.

This gives the following steady-state diagram, where the monostable values are shown as a heat map on $\zeta_{\infty}$.

Sensitivity and robustness analyses

We then turned to the bistability thresholds $\lambda_{\pm}$ defining the bistable regime. These are computed numerically using the following function.

For the following analyses, we need to generate synthetic data resulting from variability in the model parameters. Thus, we need to discuss parameter sampling. The method that we will use involves drawing parameters from log-normal distributions. Each distribution is characterized by a variation coefficient $D$, equal to the ratio of the standard deviation to the median of the distribution. Mathematically, for a specified variation $D$, the corresponding standard deviation is $$s = \sqrt{\log \left( \frac{1 + \sqrt{1 + 4 D^{2}}}{2} \right)}.$$

The cell below defines the various functions that we will need to perform the desired analyses.

We then generate $500$ parameter sets at $D = 0.25$ about baseline median values, on which we conduct the sensitivity analysis.

The $\lambda_{\pm}$ do not appear particularly robust, but there does seem to be some sense of robustness of bistability:

This suggests the following robustness analysis, where we compute the bistable fraction $f$ of parameters drawn widely over varying $D$. The data fit the exponential form $$\hat{f} = f_{\infty} + \left( 1 - f_{\infty} \right) e^{-D / D_{\mathrm{0}}},$$ where $D_{\mathrm{0}}$ is the characteristic decay length in $D$. The fitted value of $f_{\infty}$ provides evidence for a nontrivial asymptotic bistable fraction.

Cell-level cluster integration

For the cell-level integration, we draw parameters for $100$ clusters at $D = 0.25$ about baseline median values. These define heterogeneous clusters, which we integrate into a normalized measure of cell activation as $$\zeta_{\infty}^{\mathrm{cell}} = \frac{\sum_{i} \zeta_{\infty, i}}{\sum_{i} \sigma_{i}},$$ where the subscript $i$ denotes reference to cluster $i$.

The cell-level hysteresis curve is shown as follows, from which we see a graded response, despite the threshold properties of individual clusters as seen earlier. Furthermore, we plot the thresholds $\lambda_{\pm}$ on the threshold plane, which shows a significant linear dependence between the $\lambda_{\pm}$ (the valid region $\lambda_{+} > \lambda_{-}$ is shown in green).

Model discrimination

In the following, we demonstrate model discrimination between our model, hereafter termed the cluster model, and the prevailing crosslinking model using steady-state invariants. To do this, we need to discuss the crosslinking model, which we do now.

The crosslinking model that we consider has the reactions $$L + R \mathop{\longleftrightarrow}^{k_{+}}_{k_{-}} C_{1}, \quad C_{1} + R \mathop{\longleftrightarrow}^{2 k_{+}}_{2 k_{-}} C_{2}, \quad C_{2} + R \mathop{\longleftrightarrow}^{k_{+}}_{3 k_{-}} C_{3},$$ where $L$ is FasL, $R$ is Fas, and $C_{i}$ is the complex FasL:Fas$_{i}$ for $i = 1, 2, 3$. All reaction rates have been adjusted for combinatorial factors. With the nondimensionalizations $$\lambda = \frac{l}{s}, \quad \rho = \frac{r}{s}, \quad \gamma_{i} = \frac{c_{i}}{s}, \quad \kappa = \frac{k_{-}}{k_{+} s}, \quad \tau = k_{+} st,$$ the mass-action dynamics are $$\frac{d \lambda}{d \tau} = -\nu_{1}, \quad \frac{d \rho}{d \tau} = -\nu_{1} - \nu_{2} - \nu_{3}, \quad \frac{d \gamma_{1}}{d \tau} = \nu_{1} - \nu_{2}, \quad \frac{d \gamma_{2}}{d \tau} = \nu_{2} - \nu_{3}, \quad \frac{d \gamma_{3}}{d \tau} = \nu_{3},$$ where $$\nu_{1} = 3 \lambda \rho - \kappa \gamma_{1}, \quad \nu_{2} = 2 \rho \gamma_{1} - 2 \kappa \gamma_{2}, \quad \nu_{3} = \rho \gamma_{2} - 3 \kappa \gamma_{3}.$$ Therefore, at steady state, $$\gamma_{1, \infty} = 3 \lambda_{\infty} \left( \frac{\rho_{\infty}}{\kappa} \right), \quad \gamma_{2, \infty} = 3 \lambda_{\infty} \left( \frac{\rho_{\infty}}{\kappa} \right)^{2}, \quad \gamma_{3, \infty} = \lambda_{\infty} \left( \frac{\rho_{\infty}}{\kappa} \right)^{3},$$ where $$\lambda_{\infty} = \frac{\Lambda}{1 + 3 \left( \rho_{\infty} / \kappa \right) + 3 \left( \rho_{\infty} / \kappa \right)^{2} + \left( \rho_{\infty} / \kappa \right)^{3}}$$ and $\rho_{\infty}$ is given by solving $$\rho_{\infty}^{4} + \left( 3 \Lambda + 3 \kappa - \sigma \right) \rho_{\infty}^{3} + 3 \kappa \left( 2 \Lambda + \kappa - \sigma \right) \rho_{\infty}^{2} + \kappa^{2} \left( 3 \Lambda + \kappa - 3 \sigma \right) \rho_{\infty} - \kappa^{3} \sigma = 0,$$ where $$\Lambda = \lambda + \gamma_{1} + \gamma_{2} + \gamma_{3}, \quad \sigma = \rho + \gamma_{1} + 2 \gamma_{2} + 3 \gamma_{3}.$$ It is easy to show that only one nonnegative roots exists, namely, $$\rho_{\infty} = \frac{1}{2} \left[ \sqrt{\left( 3 \Lambda + \kappa - \sigma \right)^{2} + 4 \kappa \sigma} - \left( 3 \Lambda + \kappa - \sigma \right) \right].$$

As baseline parameters for the crosslinking model, we set $\sigma = 1$ and $\kappa = 0.1$. The following figure indicates that the resulting behavior is consistent with that for the cluster model, where we set $$\zeta \equiv \gamma_{1} + 2 \gamma_{2} + 3 \gamma_{3} = \sigma - \rho.$$

Hyperactive mutants

Denote hyperactive mutant Fas by $Z_{\Delta}$. To amend the cluster model for interaction with $Z_{\Delta}$, we replace the ligand-independent and -dependent cluster-stabilization reactions with $$\begin{align*} jY + kZ + \left( i - j - k \right) Z_{\Delta} &\xrightarrow{k_{s}^{\left( i \right)}} \left( j - k' \right) Y + \left( k + k' \right) Z + \left( i - j - k \right) Z_{\Delta}, \quad \begin{cases} i = 2, \dots, m,\\ j = 1, \dots, i,\\ k = 0, \dots, i - j,\\ k' = 1, \dots, j, \end{cases}\\ L + jY + kZ + \left( i - j - k \right) Z_{\Delta} &\xrightarrow{k_{l}^{\left( i \right)}} L + \left( j - k' \right) Y + \left( k + k' \right) Z + \left( i - j - k \right) Z_{\Delta}, \quad \begin{cases} i = 2, \dots, n,\\ j = 1, \dots, i,\\ k = 0, \dots, i - j,\\ k' = 1, \dots, j, \end{cases} \end{align*}$$ respectively. The corresponding velocities are now $$\begin{align*} \nu_{s} &= \sum_{i = 2}^{m} \kappa_{s}^{\left( i \right)} \sum_{j = 1}^{i} \eta^{j} \sum_{k = 0}^{i - j} \zeta^{k} \zeta_{\Delta}^{i - j - k} \sum_{k' = 1}^{j} k',\\ \nu_{l} &= \lambda \sum_{i = 2}^{n} \kappa_{l}^{\left( i \right)} \sum_{j = 1}^{i} \eta^{j} \sum_{k = 0}^{i - j} \zeta^{k} \zeta_{\Delta}^{i - j - k} \sum_{k' = 1}^{j} k'. \end{align*}$$

We compute the updated form of $d \zeta / d \tau$ below, where we set $$\zeta_{\Delta} \equiv \frac{z_{\Delta}}{s} = \Delta \overline{\sigma}$$ for $\overline{\sigma} = \sigma + \zeta_{\Delta}$ the total receptor (wildtype and mutant) concentration, and $\Delta$ the mutant population fraction.

We use Sage to perform some preliminary symbolic computations on the modified dynamics function.

The following plot shows the response curve of the active wildtype fraction $\varphi = \zeta_{\infty} / \sigma$ for various levels of $\Delta$ at $\overline{\sigma} = 1$ (stable, solid lines; unstable, dashed lines). That the response curve is sensitive to $\Delta$ under the cluster model is in contrast to that under the crosslinking model, which predicts an invariant $\varphi$-response curve, due to the absence of receptor interactions.

Steady-state invariants

The steady-state invariant for each model is defined as $d \zeta / d \tau$. However, the invariants must be expressed only in terms of experimentally accessible parameters, which we took as $(\Lambda, \sigma, \zeta_{\infty})$. Thus, for the cluster invariant, we make the identification $\lambda \mapsto \Lambda$, and for the crosslinking invariant, we identify $\rho_{\infty} \mapsto \sigma - \zeta_{\infty}$. The following cell defines the model invariants.

We implement a function for generating the experimentally accessible values $(\Lambda, \sigma, \zeta_{\infty})$ below.

As demonstration that model discrimination using steady state invariants is practical, we generate synthetic data for $\Lambda$ and $\sigma$; $\zeta_{\infty}$ is then calculated using the above function. For the cluster model, if bistability is encountered, one of the two stable values of $\zeta_{\infty}$ is chosen at random.

For each model-data pair, we fit the model invariant to the data by minimizing over the model parameter space using SLSQP. This gives the best possible fit for a given model to the data.

The results are summarized below, from which we see that the model can be correctly identified from the data.