# /// 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()