April 20th, 2017

$$ \newcommand{\esp}[1]{\mathbb{E}\left(#1\right)} \newcommand{\var}[1]{\mbox{Var}\left(#1\right)} \newcommand{\deriv}[1]{\dot{#1}(t)} \newcommand{\prob}[1]{ \mathbb{P}\!(#1)} \newcommand{\eqdef}{\mathop{=}\limits^{\mathrm{def}}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\bpsi}{\boldsymbol{\psi}} \def\pmacro{\texttt{p}} \def\like{{\cal L}} \def\llike{{\cal LL}} \def\logit{{\rm logit}} \def\probit{{\rm probit}} \def\one{{\rm 1\!I}} \def\iid{\mathop{\sim}_{\rm i.i.d.}} \def\simh0{\mathop{\sim}_{H_0}} \def\df{\texttt{df}} \def\res{e} \def\xomega{x} \newcommand{\argmin}[1]{{\rm arg}\min_{#1}} \newcommand{\argmax}[1]{{\rm arg}\max_{#1}} \newcommand{\Rset}{\mbox{$\mathbb{R}$}} \def\param{\theta} \def\setparam{\Theta} \def\xnew{x_{\rm new}} \def\fnew{f_{\rm new}} \def\ynew{y_{\rm new}} \def\nnew{n_{\rm new}} \def\enew{e_{\rm new}} \def\Xnew{X_{\rm new}} \def\hfnew{\widehat{\fnew}} \def\degree{m} \def\nbeta{d} \newcommand{\limite}[1]{\mathop{\longrightarrow}\limits_{#1}} \def\ka{k{\scriptstyle a}} \def\ska{k{\scriptscriptstyle a}} \def\kel{k{\scriptstyle e}} \def\skel{k{\scriptscriptstyle e}} \def\cl{C{\small l}} \def\Tlag{T\hspace{-0.1em}{\scriptstyle lag}} \def\sTlag{T\hspace{-0.07em}{\scriptscriptstyle lag}} \def\Tk{T\hspace{-0.1em}{\scriptstyle k0}} \def\sTk{T\hspace{-0.07em}{\scriptscriptstyle k0}} \def\thalf{t{\scriptstyle 1/2}} \newcommand{\Dphi}[1]{\partial_\pphi #1} \def\asigma{a} \def\pphi{\psi} \newcommand{\stheta}{{\theta^\star}} \newcommand{\htheta}{{\widehat{\theta}}}$$

The population approach may also be relevant for building predictive computational models of intracellular processes.

For example, consider as a case study the experiments performed by Uhlendorf *et al.*, 2012, looking at the high-osmolarity glycerol (HOG) pathway in budding yeast.

Yeast cells are exposed to osmotic shocks, i.e., sudden changes in the solute concentration of their surroundings. Signal transduction pathways, most notably the HOG pathway, provide information to the cell about the osmolarity of its environment and activate responses to deal with these stress conditions. In particular, a large set of genes is turned on and corresponding stress-responsive proteins are produced.

Let \(v(t)\) be the valve status at time \(t\) with \(v=1\) when the valve is “on” (hyperosmotic media) and \(v=0\) when it is “off” (normal media).

Then, \(v\) is a piecewise constant function, solution of \(\deriv{v}=0\) with source terms equal to 1 when the valve is turned on, and equal to -1 when it is turned off.

In the following example, the valve is “on” between 10’ and 18’ and “off” betwen 35’ and 50’:

```
valveModel <- inlineModel("
[LONGITUDINAL]
PK:
depot(target=v)
EQUATION:
ddt_v=0
")
ton <- list(amount = 1, time = c(10, 35))
toff <- list(amount = -1, time = c(18, 50))
out <- list(name = 'v', time = seq(0, 60, by=0.1))
res <- simulx(model=valveModel, output=out, treatment=list(ton,toff))
print(ggplot(res$v) + geom_line(aes(time,v)) + ylab('valve status'))
```

For modeling, we need first to take into account the lag between change in valve status and change of osmolarity \(h(t)\) of the cell environment: We assume that the medium needs one minute to change from normal to hyperosmotic when the valve is turned on and four minutes to return to normal when turned off.

The function \(h\) can then be defined by introducing inputs of value 1 at times \((t^{on}_k, 1\leq k \leq K)\) with a rate of 1, and inputs with value -1 at times \((t^{off}_k, 1\leq k \leq K)\) with rate 0.25. Here, \(t^{off}_k-t^{on}_k\) is the duration of the \(k\)th shock.

Then, the extracted nuclear Hog1 function \(s(t)\) for a given salt concentration is a solution to the following ODE: \[ \deriv{s} = \kappa\, h(t) - \gamma \, s(t). \] Lastly, the corresponding time-varying gene activation intensity \(u(t)\) is obtained by transforming the Hog1 abundance using a Hill function: \[ u(t) = \frac{(s(t) + s_0)^{n_H}}{K_d^{n_H} + (s(t) + s_0)^{n_H}}. \]

This model is implemented in `hog1_model`

.

```
hog1_model <- inlineModel("
[LONGITUDINAL]
input={kappa, gamma, nh, Kd, s0}
PK:
depot(target=h)
EQUATION:
ddt_h = 0
ddt_s = kappa*h - gamma*s
sh = (s+s0)^nh
u = sh/(Kd^nh+sh)
")
```

Assume again that the valve is “on” between 10’ and 18’ and “off” betwen 35’ and 50’. Status of the valve will be defined according to input source terms defined in lists `ton` and `toff`:

```
ton <- list(amount = 1,
rate = 1,
time = c(10, 35))
toff <- list(amount = -1,
rate = 0.25,
time = c(18, 50))
```

We will display \(h\), \(s\) and \(u\) computed between 0’ and 70’ with a given set of parameter values:

```
library("reshape2")
out <- list(name = c('h','s','u'),
time = seq(0, 70, by=0.1))
p1 <- c(kappa=0.3, gamma=0.7, nh=6, Kd=0.15, s0=0.016)
res <- simulx(model = hog1_model,
parameter = p1,
output = out,
treatment = list(ton, toff))
r <- merge(res$h,merge(res$s,res$u))
r <- melt(r, id = 'time', variable.name = 'signal')
print(ggplot(r, aes(time,value)) + geom_line(aes(colour = signal),size=1) +
ylab('activation')+ theme(legend.position=c(.9, .7)))
```

The related Hog1-induced gene expression model is given by a reaction network described by Zechner *et al.* (2012) and adapted by Gonzalez *et al.* (2013).

Let \(x_1\), \(x_2\) and \(x_4\) be the proportions of promoters that are in the *off* state, in the *on* state or bound to chromatin remodeling factors. Let \(x_3\), \(x_5\) and \(x_6\) be the concentrations of remodeling factors, messenger RNA and fluorescent yECitrine proteins.

The effective activation rate of the gene *pSTL1* is assumed to be proportional to the activation intensity \(u\) defined in the previous section with a delay \(\tau\).

Then, the system can be described by the following set of reaction rate equations:

\[ \begin{aligned} \deriv{x_1} & = c_2 x_2(t) - c_1 u(t-\tau) x_1(t) \\ \deriv{x_2} & = c_1 u(t-\tau) x_1(t) - c_2 x_2(t) + c_4 x_4(t) -c_3 x_2(t)x_3(t)\\ \deriv{x_3} & = c_4 x_4(t) -c_3 x_2(t)x_3(t)\\ \deriv{x_4} & = c_3 x_2(t)x_3(t) - c_4 x_4(t) \\ \deriv{x_5} & = c_5 x_4(t) -c_8 x_5(t)\\ \deriv{x_6} & = c_6 x_5(t) -c_7 x_6(t), \end{aligned} \] where the initial conditions at \(t=t_0\) are \(x_1(t_0)=1\), \(x_2(t_0)=x_4(t_0)=x_5(t_0)=x_6(t_0)=0\) and \(x_3(t_0)=x_{3,0}\).

The structural model combines this dynamical system with the model for the activation intensity described in the previous section. It is implemented in `hog2_model`

:

```
hog2_model <- inlineModel("
[LONGITUDINAL]
input={kappa, gamma, nh, Kd, s0, c1, c2, c3, c4, c5, c6, c7, c8, x30, tau}
PK:
depot(target=h,Tlag=tau)
EQUATION:
ddt_h = 0
ddt_s = kappa*h - gamma*s
sh = (s+s0)^nh
u = sh/(Kd^nh+sh)
x1_0 = 1
x3_0 = x30
ddt_x1 = c2*x2 - c1*u*x1
ddt_x2 = c1*u*x1 - c2*x2 + c4*x4 - c3*x2*x3
ddt_x3 = c4*x4 - c3*x2*x3
ddt_x4 = c3*x2*x3 - c4*x4
ddt_x5 = c5*x4 - c8*x5
ddt_x6 = c6*x5 - c7*x6
")
```

Let us use `simulx`

for computing \(x_6\) using the experimental design used by Gonzalez *et al.* (2012)

```
ton <- list(amount = 1,
rate = 1,
time = c(35, 65, 95, 125, 155, 185, 215, 245, 280, 341, 371, 401, 449, 479, 533, 563, 649, 683))
toff <- list(amount = -1,
rate = 0.25,
time = c(43, 73, 103, 133, 163, 193, 223, 253, 285, 349, 379, 409, 454, 486, 541, 570, 657, 691))
out <- list(name=c('x6'), time=seq(0, 800, by=0.5))
p2 <- c(c1=75, c2=2500, c3=0.00002, c4=0.04, c5=1.2, c6=800, c7=0.005, c8=0.1, x30=150, tau=10)
res <- simulx(model = hog2_model,
parameter = c(p1, p2),
output = out,
treatment = list(ton, toff))
plot(ggplot() + geom_line(data=res$x6, aes(x=time, y=x6)))
```

Gonzalez, A., Uhlendorf, J., Schaul, J., Cinquemani, E., Batt, G., and Ferrari-Trecate, G. (2013). Identification of biological models from single-cell data: a comparison between mixed-effects and momentbased inference. In Proceedings of ECC - 12th European Control Conference, pages 2652-2657.

Uhlendorf, J., Miermont, A., Delaveau, T., Charvin, G., Fages, F., Bottani, S., Batt, G., and Hersen, P. (2012). Long-term model predictive control of gene expression at the population and single-cell levels. Proceedings of the National Academy of Sciences, 109(35):14271-14276.

Zechner, C., Ruess, J., Krenn, P., Pelet, S., Peter, M., Lygeros, J., and Koeppl, H. (2012). Moment-based inference predicts bimodality in transient gene expression. Proceedings of the National Academy of Sciences, 109(21):8340-8345.