sciencehub-team/mrna-expressionPublic

# /// script
# requires-python = ">=3.12"
# dependencies = ["marimo>=0.25.1", "numpy", "matplotlib"]
# ///

import marimo

__generated_with = "0.25.1"
app = marimo.App(width="medium")


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        # How long does an mRNA dose make protein?

        An mRNA therapeutic or vaccine delivers instructions, not the protein
        itself. Once the mRNA reaches the cytoplasm, ribosomes translate it
        until it is degraded, and the protein it encodes is degraded in turn.
        The simplest model that describes single-cell transfection data has two
        first-order steps (Leonhardt et al. 2014):

        $$\frac{dm}{dt} = -\delta\, m, \qquad \frac{dp}{dt} = k\, m - \beta\, p$$

        with mRNA $m$, protein $p$, translation rate $k$, and degradation rates
        $\delta = \ln 2 / t_{1/2}^{\text{mRNA}}$ and
        $\beta = \ln 2 / t_{1/2}^{\text{protein}}$. Starting from $m_0$ at
        $t = 0$ the protein level is

        $$p(t) = \frac{k\, m_0}{\delta - \beta}\left(e^{-\beta t} - e^{-\delta t}\right).$$
        """
    )
    return


@app.cell
def _():
    import marimo as mo
    import numpy as np
    import matplotlib.pyplot as plt
    return mo, np, plt


@app.cell
def _(mo):
    mrna_half_life = mo.ui.slider(0.5, 48, step=0.5, value=10, label="mRNA half-life (h)")
    protein_half_life = mo.ui.slider(1, 96, step=1, value=26, label="Protein half-life (h)")
    mo.vstack([mrna_half_life, protein_half_life])
    return mrna_half_life, protein_half_life


@app.cell
def _(np):
    def protein(t, mrna_half, protein_half, k=1.0, m0=1.0):
        """Protein level after one dose at t = 0, in units of k * m0 * hours."""
        delta, beta = np.log(2) / mrna_half, np.log(2) / protein_half
        t = np.asarray(t, dtype=float)
        if np.isclose(delta, beta):
            level = k * m0 * t * np.exp(-beta * t)
        else:
            level = k * m0 / (delta - beta) * (np.exp(-beta * t) - np.exp(-delta * t))
        return np.where(t >= 0, level, 0.0)

    def peak_time(mrna_half, protein_half):
        delta, beta = np.log(2) / mrna_half, np.log(2) / protein_half
        if np.isclose(delta, beta):
            return 1 / beta
        return np.log(delta / beta) / (delta - beta)

    return peak_time, protein


@app.cell
def _(mo, mrna_half_life, np, peak_time, plt, protein, protein_half_life):
    hours = np.linspace(0, 168, 800)
    level = protein(hours, mrna_half_life.value, protein_half_life.value)
    mrna = np.exp(-np.log(2) / mrna_half_life.value * hours)
    t_peak = peak_time(mrna_half_life.value, protein_half_life.value)
    exposure = mrna_half_life.value * protein_half_life.value / np.log(2) ** 2

    fig, ax = plt.subplots(figsize=(7, 3.6))
    ax.plot(hours / 24, level / level.max(), color="#0969da", lw=2, label="protein")
    ax.plot(hours / 24, mrna, color="#9a6700", lw=1.5, ls="--", label="mRNA")
    ax.axvline(t_peak / 24, color="#57606a", lw=0.8)
    ax.set(xlabel="days after the dose", ylabel="fraction of peak", xlim=(0, 7), ylim=(0, 1.05))
    ax.legend(frameon=False)
    fig.tight_layout()

    mo.vstack(
        [
            fig,
            mo.md(
                f"Protein peaks **{t_peak:.1f} hours** after the dose. Total protein made "
                f"over the whole course is $k m_0 / (\\delta\\beta)$ = **{exposure:,.0f}** "
                "in units of $k m_0$ hours squared."
            ),
        ]
    )
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ## Why mRNA stability matters

        The area under the protein curve is $k m_0 / (\delta \beta)$: it is
        proportional to the mRNA half-life times the protein half-life. Doubling
        how long the mRNA survives doubles the total protein made from the same
        dose, without changing the dose. Much of the engineering of therapeutic
        mRNA (modified nucleosides such as N1-methylpseudouridine, optimised
        untranslated regions, the cap structure) works on $k$ and $\delta$.
        The curves below keep the protein half-life you chose and vary only the
        mRNA half-life.
        """
    )
    return


@app.cell
def _(np, plt, protein, protein_half_life):
    days = np.linspace(0, 10, 800)
    fig2, ax2 = plt.subplots(figsize=(7, 3.4))
    for _half, _shade in zip((2, 6, 12, 24), ("#d0d7de", "#8c959f", "#57606a", "#0969da")):
        ax2.plot(
            days,
            protein(days * 24, _half, protein_half_life.value),
            color=_shade,
            lw=2,
            label=f"mRNA half-life {_half} h",
        )
    ax2.set(xlabel="days after the dose", ylabel="protein (k m₀ h)")
    ax2.legend(frameon=False)
    fig2.tight_layout()
    fig2
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ## Repeated doses

        The model is linear, so repeated doses add up. The protein level after
        a schedule is the sum of the single-dose curves, each shifted to its
        dose time. Choose an interval and a number of doses to see whether
        expression stays above a threshold between doses.
        """
    )
    return


@app.cell
def _(mo):
    interval = mo.ui.slider(6, 96, step=6, value=24, label="Dose every (h)")
    doses = mo.ui.slider(1, 10, step=1, value=5, label="Number of doses")
    threshold = mo.ui.slider(0.05, 1.0, step=0.05, value=0.5, label="Threshold (fraction of one-dose peak)")
    mo.vstack([mo.hstack([interval, doses], justify="start"), threshold])
    return doses, interval, threshold


@app.cell
def _(
    doses,
    interval,
    mo,
    mrna_half_life,
    np,
    plt,
    protein,
    protein_half_life,
    threshold,
):
    span = interval.value * doses.value + 96
    t = np.linspace(0, span, 1600)
    single_peak = protein(np.linspace(0, 400, 4000), mrna_half_life.value, protein_half_life.value).max()
    total = sum(
        protein(t - n * interval.value, mrna_half_life.value, protein_half_life.value)
        for n in range(doses.value)
    ) / single_peak
    covered = t[(t > 0) & (t < interval.value * doses.value)]
    above = np.mean(np.interp(covered, t, total) >= threshold.value) if covered.size else 0.0

    fig3, ax3 = plt.subplots(figsize=(7, 3.4))
    ax3.plot(t / 24, total, color="#0969da", lw=2)
    ax3.axhline(threshold.value, color="#cf222e", lw=0.8, ls="--")
    for _n in range(doses.value):
        ax3.axvline(_n * interval.value / 24, color="#d0d7de", lw=0.8)
    ax3.set(xlabel="days", ylabel="protein / one-dose peak")
    fig3.tight_layout()
    mo.vstack(
        [
            fig3,
            mo.md(
                f"During dosing, expression is above the threshold **{above:.0%}** of the time. "
                f"The highest level reached is **{total.max():.2f}** times a single dose's peak."
            ),
        ]
    )
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ---
        **Model and limits.** One compartment, first-order kinetics, all mRNA
        released at $t = 0$. Real delivery is spread over hours, cells differ
        widely in all three rates, and the immune response to a vaccine is not
        modelled. Levels are relative; set $k$ and $m_0$ from your own data to
        get absolute amounts. The default protein half-life, 26 hours, is the
        value reported for EGFP in mammalian cells (Corish and Tyler-Smith 1999).
        **References.** A. Leonhardt et al., *Nanomedicine: NBM* 10 (2014) 679;
        P. Corish and C. Tyler-Smith, *Protein Eng.* 12 (1999) 1035;
        K. Karikó et al., *Mol. Ther.* 16 (2008) 1833.
        Fork it to fit the model to your own time-course data.
        """
    )
    return


if __name__ == "__main__":
    app.run()