introduction

A few months ago I blogged a numerical simulation of prion protein degradation. Now, at Sina Ghaemmaghami‘s suggestion, I’m returning to this issue to figure out an analytical solution.

Here’s a bit more background for math people who aren’t biologists reading this post. Prion diseases are caused by a normal protein called cellular prion protein (PrPC) misfolding into a diseased conformation called scrapie prion protein (PrPSc).  In some anti-prion drug discovery screens, people take prion-infected cells and screen a whole library of tens of thousands of chemical compounds to look for compounds that will reduce the amount of PrPSc the cells accumulate.  If your library has a compound that directly reduces PrPSc, you’re pretty likely to see that effect.

But what if there’s a compound whose mechanism is way far upstream? Remember that DNA → RNA → protein and (in this case), healthy protein → diseased protein. So for instance, what if a compound inhibits a transcription factor, preventing the prion protein gene’s DNA from being transcribed into RNA?  In this case, the pre-existing messenger RNA (mRNA) will first need to decay with a certain half-life, and then as the level of mRNA gradually reduces, the amount of PrPC translated form that mRNA will gradually reduce, and then finally the amount of PrPSc converted from PrPC will gradually reduce. So switching off the supply of mRNA could take days to percolate through the system and show you an appreciable reduction in PrPSc. How many days? That’s what this post seeks to answer.

acknowledgements

When I set out to do this, it had been over a decade since I’d solved a differential equation. Thanks to Adam Bliss for a lot of help figuring out how to go about doing all this.

assumptions

I’ll assume that the degradation of PrP mRNA, PrPC and PrPSc is first-order, so they follow the standard rules of exponential decay. This is also what others have assumed and it’s pretty well-supported by experimental evidence [Ghaemmaghami 2007 + supplement].

The decay parameters of mRNA, PrPC and PrPSc will be denoted respectively by λ, μ and ν. I’ll assume the production rate of PrP mRNA is a constant denoted r. The production rate of PrPC is a constant a times the amount of mRNA present. The production rate of PrPSc is a constant β times the amount of PrPC present*. The quantities of PrP mRNA, PrPC and PrPSc as a function of time will be denoted by R(t), S(t) and T(t) respectively. The rates of change in these variables are then described in a system of differential equations:

\text{R'(t)} = r - \lambda\text{R(t)}
\text{S'(t)} = \alpha\text{R(t)} - \mu\text{S(t)}
\text{T'(t)} = \beta\text{S(t)} - \nu\text{T(t)}

*I’ve assumed that the rate of PrPSc production is proportional to the amount of PrPC available as substrate but does not depend on the amount of PrPSc available to do the converting. This is surely a controversial assumption, as prion conversion requires PrPSc-PrPC binding. I’ll return to this subject in detail in an upcoming post. For now, I will stick to the assumption that the PrPSc formation rate is independent of [PrPSc], as I assumed in my original numerical simulations.

Finally, I’ll also assume we begin at steady state at t = 0 and will define the initial steady state quantities as:

R_0 = R(0)
S_0 = S(0)
T_0 = T(0)

For simplicity, when I plot later I’ll assume that all of those quantities are equal to 1.

Based on the literature reviewed in my original post the half-lives of mRNA, PrPC and PrPSc are about 7h, 5h and 30h, so I will define the constants λ0, μ0, and ν0 to be ln(2)/7 = .099h-1, ln(2)/5 = .14h-1, and ln(2)/30 = .023h-1 respectively. I will consider these as constants while I will use the symbols λ, μ, and ν to denote variables which can be modified by chemical compounds to perturb the cell from its initial steady state.

derivation

In my previous post I considered six possible mechanisms by which small molecules could reduce PrPSc accumulation. They might inhibit transcription (reducing r), promote mRNA degradation (increasing λ), inhibiting translation (reducing α), promoting PrPC degradation (increasing μ), inhibiting conversion (reducing β) or promoting PrPSc degradation (increasing ν).

First, let’s characterize the initial steady state at t = 0, before the cell has been perturbed by any chemicals. Since I’m going to treat r, α and β as variables that can be perturbed, I’ll need to define the constants r0, α0 and β0 to be the initial, unperturbed values of these variables. At t = 0,

R'(0) = 0 = r_0 - \lambda_0 R(0)
R_0 = \frac{r_0}{\lambda_0}
S'(0) = 0 = \alpha_0 R(0) - \mu_0 S(0)
S_0 = \frac{\alpha_0 R_0}{\mu_0} = \frac{\alpha_0 r_0}{\mu_0 \lambda_0}
T'(0) = 0 = \beta_0 S(0) - \nu_0 T(0)
T_0 = \beta_0 S_0 / \nu_0 = \frac{\beta_0 \alpha_0 r_0}{\nu_0 \mu_0 \lambda_0}

From these equations you can also probably guess a more general fact about the steady state values of all of these variables. For example, if I cut transcription by ½ or if I double the rate of mRNA degradation then I’ll wind up with half as much mRNA, and therefore half as much PrPC and half as much PrPSc as before.

My goal here is to see how quickly the system will advance toward that new steady state if we start from the initial steady state I’ve defined above. In other words, I need to obtain expressions for R(t), S(t) and T(t). Let’s start with R(t). I can get that by solving this differential equation:

R'(t) = r - \lambda R(t)

Which I’ll do by guess-and-check. I’ll guess that R(t) has the form A + Be-λt. Let’s try it out:

R(t) = A + Be^{-\lambda t}
R'(t) = r - \lambda (A + Be^{-\lambda t})
-\lambda Be^{-\lambda t} = r - \lambda (A + Be^{-\lambda t})
0 = r - \lambda A + \lambda Be^{-\lambda t} - \lambda Be^{-\lambda t}
A = \frac{r}{\lambda}

By gathering the constant terms together we can see that A = r/λ. Unfortunately we can’t learn anything by gathering the e-λt terms because they cancel. However we can take advantage of the assumption that we are at steady state at t = 0 to observe that:

R(0) = R_0 = A + Be^{0}
R_0 = A + B
B = R_0 - A
B = R_0 - \frac{r}{\lambda}

So now we can define B in terms of our starting point, R0. Plugging A and B back into the equation for R(t) we get:

R(t) = \frac{r}{\lambda} + (R_0 - \frac{r}{\lambda})e^{-\lambda t}

This makes perfect sense – this is just a classic exponential decay situation with one twist. r/λ represents the new steady state that we’ll asymptotically approach – think of it as the mRNA that is being replaced. (R0 – r/λ) represents the portion of the initial mRNA that is not being replaced – that portion will decay exponentially with parameter λ

And let’s just double check that it gives us the derivative we were hoping for:

\frac{dR}{dt} = -\lambda (R_0 - \frac{r}{\lambda})e^{-\lambda t}

That looks pretty different from the original expression we had for the derivative:

\frac{dR}{dt} = r - \lambda R(t)

But in fact they’re the same:

\frac{dR}{dt} = r - \lambda (\frac{r}{\lambda} + (R_0 - \frac{r}{\lambda})e^{-\lambda t})
\frac{dR}{dt} = r - \lambda \frac{r}{\lambda} - \lambda(R_0 - \frac{r}{\lambda})e^{-\lambda t})
\frac{dR}{dt} = - \lambda(R_0 - \frac{r}{\lambda})e^{-\lambda t})

Moving on we can apply the same bag of tricks to get S(t). Let’s guess that S(t) has the form A + Be-λt + Ce-μt.

S(t) = A + Be^{-\lambda t} + Ce^{-\mu t}
S'(t) = \alpha R(t) - \mu S(t)
S'(t) = -\lambda Be^{-\lambda t} - \mu Ce^{-\mu t}
-\lambda Be^{-\lambda t} - \mu Ce^{-\mu t} = \alpha (\frac{r}{\lambda} + (R_0 - \frac{r}{\lambda})e^{-\lambda t}) - \mu (A + Be^{-\lambda t} + Ce^{-\mu t})
0 = \alpha \frac{r}{\lambda} - \mu A + \lambda Be^{-\lambda t} + (R_0 - \frac{r}{\lambda})e^{-\lambda t}) - \mu Be^{-\lambda t} + \mu Ce^{-\mu t} - \mu Ce^{-\mu t}

Gathering the like terms we get:

A = \frac{\alpha r}{\mu \lambda}
B = \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}

And this time, we can’t learn anything about C except by observing that at t = 0,

S(0) = S_0 = A + B + C
C = S_0 - A - B
C = S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}

Putting it all together we find that:

S(t) = \frac{\alpha r}{\mu \lambda} + (\alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\lambda t} + (S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\mu t}

And similarly we’ll guess that T(t) has the form A + Be-λt + Ce-μt + De-νt.

T(t) = A + Be^{-\lambda t} + Ce^{-\mu t} + De^{-\nu t}
T'(t) = \beta S(t) - \nu T(t)
T'(t) = -\lambda Be^{-\lambda t} - \mu Ce^{-\mu t} - \nu De^{-\nu t}
T'(t) = -\lambda Be^{-\lambda t} - \mu Ce^{-\mu t} - \nu De^{-\nu t} = \beta (\frac{\alpha r}{\mu \lambda} + (\alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\lambda t} + (S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\mu t}) - \nu (A + Be^{-\lambda t} + Ce^{-\mu t} + De^{-\nu t})

And rearranging to group together the constant, e-λt, e-μt, and e-νt terms we get:

A = \frac{\beta \alpha r}{\nu \mu \lambda}
B = \beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)}
C = \beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu}

And solving for D:
T_0 = A + B + C + D
D = T_0 - A - B - C
D = T_0 - \frac{\beta \alpha r}{\nu \mu \lambda} - \beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)} - \beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu}

Putting them all together,

T(t) = \frac{\beta \alpha r}{\nu \mu \lambda} + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)})e^{-\lambda t} + (\beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})e^{-\mu t} + (T_0 - \frac{\beta \alpha r}{\nu \mu \lambda} - \beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)} - \beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})e^{-\nu t}

It’s a good thing there’s nothing past T(t) because this is already overflowing my 600px of width on this blog. Just to review where we’re at:

R(t) = \frac{r}{\lambda} + (R_0 - \frac{r}{\lambda})e^{-\lambda t}
S(t) = \frac{\alpha r}{\mu \lambda} + (\alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\lambda t} + (S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\mu t}
T(t) = \frac{\beta \alpha r}{\nu \mu \lambda} + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)})e^{-\lambda t} + (\beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})e^{-\mu t} + (T_0 - \frac{\beta \alpha r}{\nu \mu \lambda} - \beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)} - \beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})e^{-\nu t}

Now, the one other thing we need to worry about is that my expressions for S(t) have μ – λ, ν – μ, and ν – λ in denominators. That means that if any two of λ, μ and ν are equal, then we’re dividing by zero. Intuitively, there’s nothing wrong with having the same half-life for, say, PrP mRNA and PrPC, so we should be able to find solutions for all of these cases.

First, let’s consider S(t) in the case where λ = μ. I’ll rearrange our expression for S(t) above:

S(t) = \frac{\alpha r}{\mu \lambda} + (\alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\lambda t} + (S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\mu t}
S(t) = \frac{\alpha r}{\mu \lambda} + (\alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) (e^{-\lambda t} - e^{-\mu t}) + (S_0 - \frac{\alpha r}{\mu \lambda}) e^{-\mu t}
S(t) = \frac{\alpha r}{\mu \lambda} + \alpha (R_0 - \frac{r}{\lambda})\frac{e^{-\lambda t} - e^{-\mu t}}{\mu - \lambda} + (S_0 - \frac{\alpha r}{\mu \lambda}) e^{-\mu t}

Now I have all of the μ and λ “issues” concentrated in the second term. To apply l’Hopital’s rule I need to treat t as a constant, and μ as a variable, and take the limit as μ approaches λ of derivative of this term with respect to μ:

\lim_{\mu \to \lambda} \frac{d}{d\mu} \alpha ( R_0 - \frac{r}{\lambda})\frac{e^{-\lambda t} - e^{-\mu t}}{\mu - \lambda}
\lim_{\mu \to \lambda} \alpha (R_0 - \frac{r}{\lambda}) --t \frac{e^{-\mu t}}{1}
\lim_{\mu \to \lambda} \alpha (R_0 - \frac{r}{\lambda})t e^{-\mu t}
\alpha ( R_0 - \frac{r}{\mu})t e^{-\mu t}

Now plugging this term back into the whole equation and setting λ = μ everywhere, we get:

S(t) = \frac{\alpha r}{\mu^2} + (\alpha R_0 - \frac{r}{\mu})t e^{-\mu t} + (S_0 - \frac{\alpha r}{\mu^2}) e^{-\mu t}

T(t) is an even bigger pain, since we have four possible combinations: λ = μ, μ = ν, ν = λ, and λ = μ = ν. I’ll start with ν = μ which is exactly analogous to what we just did. It will help if we rearrange the expression for T(t) as follows:

T(t) = T_0 e^{-\nu t} + \frac{\beta \alpha r}{\nu \mu \lambda}(1 - e^{-\nu t}) + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)})(e^{-\lambda t} - e^{-\nu t}) + (\beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})(e^{-\mu t} - e^{-\nu t})

Let’s consider that last term as follows:

\beta(S_0 - \frac{\alpha r}{\mu \lambda} - \alpha\frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda})(\frac{e^{-\mu t} - e^{-\nu t}}{\nu - \mu})

Taking the limit of the derivative:

\lim_{\nu \to \mu} \beta(S_0 - \frac{\alpha r}{\mu \lambda} - \alpha\frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) \frac{d}{d\nu} (\frac{e^{-\mu t} - e^{-\nu t}}{\nu - \mu})
\lim_{\nu \to \mu} \beta(S_0 - \frac{\alpha r}{\mu \lambda} - \alpha\frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) te^{-\nu t}

Plugging this back into the overall equation and setting μ to ν everywhere:

T(t) = T_0 e^{-\nu t} + \frac{\beta \alpha r}{\nu^2 \lambda}(1 - e^{-\nu t}) + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)^2})(e^{-\lambda t} - e^{-\nu t}) + \beta(S_0 - \frac{\alpha r}{\nu \lambda} - \alpha\frac{R_0 - \frac{r}{\lambda}}{\nu - \lambda}) te^{-\nu t}

If we want to go a step further and consider when λ = μ = ν, then we need to take that expression and take the derivative with respect to ν as ν → λ. Multiplying through:

T(t) = T_0 e^{-\nu t} + \frac{\beta \alpha r}{\nu^2 \lambda}(1 - e^{-\nu t}) + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)^2})(e^{-\lambda t} - e^{-\nu t}) + \beta S_0 te^{-\nu t} - \beta\frac{\alpha r}{\nu \lambda}te^{-\nu t} - \beta\alpha\frac{R_0 - \frac{r}{\lambda}}{\nu - \lambda} te^{-\nu t}

We have just two terms which pose problems. Applying l’Hopital’s rule to the final term we get:

\lim_{\nu \to \lambda} \frac{d}{d\nu} \beta \alpha (R_0 - \frac{r}{\lambda})\frac{te^{-\nu t}}{\nu - \lambda}

The derivative of the denominator is just 1. The derivative of the numerator is, by the product rule, t(-t)e-νt + (1)e-νt. Thus this term becomes:

\beta \alpha (R_0 - \frac{r}{\nu}) (1 -t^2)e^{-\nu t}

Our other problem is the third term in T(t), which contains (ν – λ)2 in the denominator. Applying l’Hopital’s rule once we can convert that fraction from:

\frac{e^{-\lambda t} - e^{-\nu t}}{(\nu - \lambda)^2}

To:

\frac{te^{-\nu t}}{2\nu - 2\lambda}

Which is still indeterminate when ν = λ. Thus we need to apply the rule a second time, as shown in the Wikipedia examples, and get:

\frac{(1-t^2)e^{-\nu t}}{2}

Conveniently, this looks similar to the other term we just applied the rule to. Plugging these terms back in and setting λ = ν everywhere,

T(t) = T_0 e^{-\nu t} + \frac{\beta \alpha r}{\nu^3}(1 - e^{-\nu t}) + \beta \alpha (R_0 - \frac{r}{\nu})\frac{(1-t^2)e^{-\nu t}}{2} + \beta S_0 te^{-\nu t} - \beta\frac{\alpha r}{\nu^2}te^{-\nu t} - \beta\alpha(R_0 - \frac{r}{\nu})(1 -t^2)e^{-\nu t}
T(t) = \frac{\beta\alpha r}{\nu^3} + e^{-\nu t}(T_0 - \frac{\beta\alpha r}{v^3} - \frac{1}{2}\beta\alpha(R_0 - \frac{r}{\nu})(1-t^2) + \beta(S_0 - \frac{\alpha r}{v^2})t)

The case where ν = λ but ν ≠ μ can be solved similarly to the ν = μ case. The case where μ = λ is even more complicated and requires you go back to the special version of S(t) for when μ = λ and re-solve for T(t). It all gets unbelievably tedious and I won’t get into it here.

Summarizing what we have so far:

R(t) = \frac{r}{\lambda} + (R_0 - \frac{r}{\lambda})e^{-\lambda t}
<br /> S(t) = \left\{<br /> \begin{array}{ll}<br /> \frac{\alpha r}{\mu^2} + \alpha(R_0 - \frac{r}{\mu})t e^{-\mu t} + (S_0 - \frac{\alpha r}{\mu^2}) e^{-\mu t} & \mbox{if } \mu = \lambda \\<br /> \frac{\alpha r}{\mu \lambda} + (\alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\lambda t} + (S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}) e^{-\mu t} & \mbox{ otherwise } \end{array}<br /> \right.<br />
<br /> T(t) = \left\{<br /> \begin{array}{ll}<br /> \frac{\beta \alpha r}{\nu \mu \lambda} + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)})e^{-\lambda t} + (\beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})e^{-\mu t} + (T_0 - \frac{\beta \alpha r}{\nu \mu \lambda} - \beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)(\mu - \lambda)} - \beta \frac{S_0 - \frac{\alpha r}{\mu \lambda} - \alpha \frac{R_0 - \frac{r}{\lambda}}{\mu - \lambda}}{\nu - \mu})e^{-\nu t} & \mbox{if } \mu \neq \lambda \;\&\; \mu \neq \nu \;\&\; \nu \neq \lambda \\<br /> \frac{\beta\alpha r}{\nu^3} + e^{-\nu t}(T_0 - \frac{\beta\alpha r}{v^3} - \frac{1}{2}\beta\alpha(R_0 - \frac{r}{\nu})(1-t^2) + \beta(S_0 - \frac{\alpha r}{v^2})t) & \mbox{if } \lambda = \mu = \nu \\<br /> T_0 e^{-\nu t} + \frac{\beta \alpha r}{\nu^2 \lambda}(1 - e^{-\nu t}) + (\beta \alpha \frac{R_0 - \frac{r}{\lambda}}{(\nu - \lambda)^2})(e^{-\lambda t} - e^{-\nu t}) + \beta(S_0 - \frac{\alpha r}{\nu \lambda} - \alpha\frac{R_0 - \frac{r}{\lambda}}{\nu - \lambda}) te^{-\nu t} & \mbox{if } \mu = \nu \;\&\; \mu \neq \lambda \\<br /> yet\;other\;stuff & \mbox {otherwise}<br /> \end{array}<br /> \right.<br />

Now that I have my analytical solution, I want to check that it matches the simulation I created a few months ago. I put the simulated and analytical solutions side by side in this R script:

# Eric Minikel
# CureFFI.org
# 2013-12-03
# Comparison of numerically simulated and analytically derived PrP degradation models

# half lives from literature
mrna_thalf  =  7 # Pfeiffer 1993  http://www.ncbi.nlm.nih.gov/pubmed/8095862
prpc_thalf  =  5 # Borcheldt 1990 http://www.ncbi.nlm.nih.gov/pubmed/1968466/
prpsc_thalf = 30 # Peretz 2001    http://www.ncbi.nlm.nih.gov/pubmed/11507642
division_time = 24 # Ghaemmaghami 2007  http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2084281/

# values at initial steady state
l0 = log(2) / mrna_thalf  # lambda_0
m0 = log(2) / prpc_thalf  # mu_0
v0 = log(2) / prpsc_thalf # nu_0
xi = log(2) / division_time

hrs = 8*24 # number of hours to model and plot. here, set to 8 days

# values of parameters, which you can tweak to explore the effects of diff. compounds
r = (l0+xi) * 0.50 # try setting to 0.50 to reduce transcription rate by half
l = (l0+xi) * 1.00 # try setting to 2.00 to double the mRNA degradation rate
a = (m0+xi) * 1.00
m = (m0+xi) * 1.00 
b = (v0+xi) * 1.00
v = (v0+xi) * 1.00

# initial steady state
R0 = 1
S0 = 1
T0 = 1

########### simulation #################

mrna_produced = numeric(hrs)
mrna_degraded = numeric(hrs)
mrna_present = numeric(hrs)
mrna_present[1] = R0
mrna_produced = rep(r, hrs)
for (i in 2:hrs) {
   mrna_degraded[i] = l * mrna_present[i-1]
   mrna_present[i] = mrna_present[i-1] + mrna_produced[i] - mrna_degraded[i]
}

prpc_produced = numeric(hrs)
prpc_degraded = numeric(hrs)
prpc_present = numeric(hrs)
prpc_present[1] = S0
prpc_produced = a * mrna_present
for (i in 2:hrs) {
   prpc_degraded[i] = m * prpc_present[i-1]
   prpc_present[i] = prpc_present[i-1] + prpc_produced[i] - prpc_degraded[i]
}

prpsc_produced = numeric(hrs)
prpsc_degraded = numeric(hrs)
prpsc_present = numeric(hrs)
prpsc_present[1] = T0
prpsc_produced = b * prpc_present
for (i in 2:hrs) {
   prpsc_degraded[i] = v * prpsc_present[i-1]
   prpsc_present[i] = prpsc_present[i-1] + prpsc_produced[i] - prpsc_degraded[i]
}

plottitle = 'Degradation of PrP mRNA, PrPC and PrPSc over time \n 50% transcription reduction, 30% translation reduction and 2x PrPSc degradation'

png('degradation.model.png',width=600,height=400)

# plot mRNA, PrPC and PrPSc values from simulation
plot(1:hrs,mrna_present,pch=1,col='red',ylim=c(0,1),ylab='normalized level',xlab='hours of incubation',main=plottitle,cex.main=.8,xaxt='n')
points(1:hrs,prpc_present,pch=1,col='orange')
points(1:hrs,prpsc_present,pch=1,col='purple')

######### analytical solutions #############

# to compare to simulated solutions, plot 1:hrs vs. 0:(hrs-1) because
# in simulation, knockdown begins at t = 1 hr, while in analytical
# solution, knockdown begins at t = 0 hr

mrna = function(t) {
  return ( r/l + (R0 - r/l)*exp(-l*t) )
}
points(1:hrs, mrna(0:(hrs-1)), type='l', lwd=2, col='red')

prpc = function(t) {
  if (l == m) {
    return ( a*r/(m^2) + a*(R0 - r/m)*t*exp(-m*t) + (S0 - a*r/(m^2))*exp(-m*t) )
  } else {
    return ( a*r/(m*l) + a*((R0 - r/l)/(m-l))*exp(-l*t) + (S0 - a*r/(m*l) - a*(R0 - r/l)/(m-l))*exp(-m*t) )
  }
}
points(1:hrs, prpc(0:(hrs-1)), type='l', lwd=2, col='orange')

prpsc = function(t) {
  if (l == m && m == v) {
    return ( b*a*r/(v^3) + exp(-v*t)*(T0 - b*a*r/(v^3) - b*a*(R0-r/v)*(1-t^2)/2 + b*(S0-a*r/(v^2))*t ) )
  } else  if (m == v) {
    return ( T0*exp(-v*t) + (b*a*r/(l*v^2))*(1-exp(-v*t)) + (b*a*(R0-r/l)/((v-l)^2))*(exp(-l*t)-exp(-v*t)) +
             b*(S0 - a*r/(v*l) - a*(R0 - r/l)/(v-l))*t*exp(-v*t) )
  } else if (l == m) {
    stop("I didn't solve those equations.")
  } else if (l == v) {
    stop("I didn't solve those equations.")
  } else {
    return ( b*a*r/(v*m*l) + b*a*(R0-r/l)/((v-l)*(m-l))*exp(-l*t) + b*(S0 - a*r/(m*l) - a*(R0-r/l)/(m-l))/(v-m)*exp(-m*t) + 
           (T0 - b*a*r/(v*m*l) - b*a*(R0-r/l)/((v-l)*(m-l)) - b*(S0 - a*r/(m*l) - a*(R0-r/l)/(m-l))/(v-m))*exp(-v*t) )
  }
}
points(1:hrs, prpsc(0:(hrs-1)), type='l', lwd=2, col='purple')

legend('topright',c('PrPSc simulated','PrPC simulated','mRNA simulated','PrPSc analytical','PrPC analytical','mRNA analytical'),
       col=c('purple','orange','red','purple','orange','red'),
       pch=c(1,1,1,NA,NA,NA), lwd=c(NA,NA,NA,2,2,2))

dev.off()

And the curves match up almost perfectly:

The very slight difference between the circles (simulation) and lines (analytical) is, I believe, analogous to the difference between “period compounding” (simulation) and “continuous compounding” (analytical). If that’s true, then making the compounding periods very fine-grained in the simulation should make the lines closer together. Rather than re-write the whole script to re-define the intervals as 15 minutes instead of 1 hour each, I just tried plugging in new half lives that are 4 times longer, since this is equivalent. And indeed, the simulation and analytical solution are even tighter in this plot:

In the above plots I’ve only considered tweaking one variable – for instance, reducing r with a compound that inhibits PrP mRNA transcription. In chemical screens in real life we’ll only test one compound at a time, but if my models (both simulated and analytical) are correct, they should also be able to predict the kinetics of treating cells with a drug cocktail that, say reduces transcription by half, reduces translation by 30% and doubles the PrPSc degradation rate. Just to check that I’ve got my equations right, let’s try out that scenario:

Once again, the simulation and analytical solution agree perfectly. And notice that mRNA’s new steady state is 50% of its original level, PrPC‘s is .5*.7 = 35% of its original level, and PrPSc‘s is .5 * .7 / 2 = 17.5% of its original level.

interpretation

Now, my original purpose in doing the simulation and now in solving these equations was to determine how long you need to incubate cells with candidate antiprion compounds in order to be able to detect compounds with a particular mechanism of action. In published high-throughput screening efforts, people have checked PrPC after 1 day, PrPSc after 5 days or PrPSc after 6 days. Is that enough? I explored this question with plots in my earlier post, but now we can ask the question a bit more rigorously.

We can plug in t1 = 1*24, 5*24 or 6*24 into the equations for R(t), S(t) and T(t) and compare the values of those functions at t1 versus at t = ∞. 

########### what percent of possible degradation is achieved by 1 day, 5 days and 6 days?

maxint = .Machine$integer.max # longest possible time in future to model

# transcription inhibition
r = l0 * 0.50 # 
l = l0 * 1.00 # 
a = m0 * 1.00
m = m0 * 1.00 
b = v0 * 1.00
v = v0 * 1.00
# percent of maximum inhibition achived after various time periods
(1 - prpc(1*24))/(1-prpc(maxint))   # .76
(1 - prpsc(5*24))/(1-prpsc(maxint)) # .90
(1 - prpsc(6*24))/(1-prpsc(maxint)) # .94

We see that, for instance, for compounds inhibiting PrP transcription, only 76% of maximal inhibition in PrPC is achieved by 24 hours.

Though it’s a bit difficult to see from examining the equations for R(t), S(t) and T(t), by playing with the above piece of code I also figured out that these half maximal inhibition values for reducing r, α or β are independent of how much you knock down r, α or β. So if you reduce r by 50%, then you get a 50% reduction in mRNA at t = ∞ and a 50%*76% = 38% reduction at t = 24 hours, just as if you reduce r by 100% (completely shutting off transcription) you get a 100% reduction at t = ∞ and a 76% reduction by t = 24 hours.

For λ, μ and ν, it’s not so simple. If you double λ you will eventually get to a steady state with 50% the original amount of mRNA, and by 24 hours you’ll be at 90% of maximal inhibition. If you multiply λ by 10, then not only will you reach a lower steady state (10% of you original mRNA level), you’ll also get there faster, with 96% of maximal inhibition by t = 24 hours.

discussion

In this post I set out to derive equations to match the simulation I had obtained earlier, and it looks like I’ve achieved that. This still leaves plenty of room open for debate as to whether I have chosen the right model to begin with. The controversial points will be that I haven’t modeled cell division, and I have considered the rate of PrPSc formation to be independent of the amount of PrPSc. I’ll plan to discuss these issues more in a future post.