The Laplace Transform in Biology

Arthur Mattuck opened his Laplace transform lectures in MIT’s 18.03 by admitting that most students meet the transform as a machine with no visible origin. A short clip of that opening, Where the Laplace Transform comes from, is the prompt for this essay. The clip covers the first minutes of the argument. The full lecture is Lecture 19 of 18.03 on MIT OpenCourseWare, and the second half of the derivation below follows the transcript published there. This piece is part of a series for quantitative readers on how molecular biology gets turned into mathematics. The first half explains where the transform comes from. The second half shows what it does in five fields of biology. It closes with what a linear-systems view means for a person’s own glucose and blood-marker time series.

A power series with a continuous index

Mattuck starts from something every reader has seen. A power series is a sum $A(x) = a(0) + a(1)\,x + a(2)\,x^2 + \cdots$, or $A(x) = \sum_{n=0}^{\infty} a(n)\,x^n$ with $n$ running from zero to infinity. He writes the coefficients as $a(n)$ rather than the usual subscript on purpose. The coefficients are a function of the discrete variable $n$: they assign a real number to each nonnegative integer. Summing the series is then an operation that takes a discrete function $a(n)$ and returns a continuous function $A(x)$.

He gives two examples. If $a(n) = 1$ for every $n$, the series is $1 + x + x^2 + x^3 + \cdots$, which he tells the class “you are supposed to be born knowing” adds up to $\frac{1}{1 - x}$. The catch is that the sum only exists when $x$ lies strictly between -1 and 1. If $a(n) = 1/n!$, the series is the exponential $e^x$ and it converges for every $x$. The clip ends where the point lands: what goes in is a function of $n$ and what comes out is a function of $x$. That change of variable is the whole idea.

The rest of the derivation comes from the OpenCourseWare transcript. Make the index continuous. Replace $n$ by a real variable $t$ running from zero to infinity, and replace the sum by the only thing that can play its role over the reals, an integral. The continuous analogue of the power series is $\int_0^\infty a(t)\,x^t\,dt$. For each fixed $x$ this integral is a number, so the result is again a function of $x$.

Two cosmetic changes turn this into the familiar form. Nobody wants $x$ as the base of an exponential, because $e$ is the only base anyone likes to differentiate or integrate. So write $x^t = e^{t \ln x}$. Then look at convergence. The integral has a chance of converging only when $x$ is between 0 and 1, so that large powers of $x$ shrink. In that range $\ln x$ is negative, and it is more pleasant to work with a positive parameter, so define $s = -\ln x$ with $s > 0$. Finally rename $a(t)$ as $f(t)$ and the output as $F(s)$. The result is

$$F(s) = \int_0^\infty f(t)\,e^{-st}\,dt$$

and this is the Laplace transform. In Mattuck’s words it is “nothing more than the continuous analog of the summation of a power series”. The convergence condition on $x$ has become a convergence condition on $s$, and the mysterious factor $e^{-st}$ is just $x^t$ after a change of name.

He draws out two features. First, this is a transform and not an operator. An operator such as differentiation takes a function of $t$ and returns a function of $t$. A transform changes the variable: a function of time $t$ goes in and a function of $s$ comes out, and in applications $s$ is very often a frequency. Second, the transform is linear, because the integral is linear. The first worked example is the transform of the constant function 1, which is $1/s$. The calculation only works when $s > 0$, and he insists that asking for the value at negative $s$ is a meaningless question, exactly as asking for the sum of the geometric series at $x = 2$ is meaningless.

Why the transform earns its keep

Two standard rules, which the clip does not reach, are why engineers reach for the transform. Integration by parts gives the derivative rule $\mathcal{L}[f'](s) = s\,F(s) - f(0)$: differentiating in time becomes multiplying by $s$ and subtracting the initial value. A linear differential equation with constant coefficients therefore becomes an algebraic equation in $s$. The second rule concerns convolution. If a system responds to an impulse with $h(t)$ and is driven by an input $u(t)$, its output is the convolution $y(t) = \int_0^t h(t - v)\,u(v)\,dv$, which is a sum of shifted and scaled impulse responses. In the Laplace domain that integral collapses to a product, $Y(s) = H(s)\,U(s)$. The function $H(s)$ is called the transfer function, and it characterises the system independently of any particular input. To get back to time you invert the transform, usually by partial fractions and a table.

Most quantitative models in biology are linear differential equations with constant coefficients or are linearised around a steady state to become one. Each such model has a transfer function, and several fields have built their standard methods on that fact.

Pharmacokinetics: compartments in the s-domain

Pharmacokinetics describes what the body does to a drug, and its workhorse is the compartment model. The body is treated as a small number of well-mixed pools, and transfer between pools and elimination from them are first-order: the rate is proportional to the amount present. The simplest case is one compartment with an intravenous bolus dose $D$. The amount $A(t)$ obeys $\frac{dA}{dt} = -kA$ with $A(0) = D$. Applying the derivative rule gives $s\,A(s) - D = -k\,A(s)$, so $A(s) = \frac{D}{s + k}$, and the table gives $A(t) = D\,e^{-kt}$.

The method scales. In a two-compartment model the transform of the central-compartment amount is a ratio of polynomials in $s$. Its denominator factors as $(s + \alpha)(s + \beta)$, so the inverse transform is the biexponential $C(t) = P\,e^{-\alpha t} + Q\,e^{-\beta t}$ that pharmacokineticists fit to plasma data. The input enters as a separate factor. A bolus is a constant $D$ in the transform domain, a constant infusion at rate $R$ is $R/s$, and first-order absorption from the gut with rate $k_a$ is $\frac{k_a D}{s + k_a}$. Concentration is then $C(s) = H(s)\,U(s)$: the same disposition function $H(s)$ multiplied by whichever input transform applies. Gibaldi and Perrier’s textbook devotes a chapter to the method of Laplace transforms and uses it throughout for linear models ¹. Benet gave the general solution for linear mammillary models, which are models with a central compartment joined to peripheral ones like spokes, with elimination allowed from any compartment ².

The linear assumption is also where the method stops. When elimination saturates, so that the enzyme clearing the drug approaches its maximum rate, the rate becomes $\frac{V_{\max} C}{K_m + C}$ rather than $kC$. The transform of a nonlinear function of $C$ is not a function of $C(s)$, and the algebra no longer closes. Gibaldi and Perrier treat this case in a separate chapter on nonlinear pharmacokinetics ³.

Gene circuits as transfer functions

Synthetic biology builds circuits from genes, and control theory supplies the language for what those circuits do to signals. Del Vecchio and Murray’s textbook develops that language for biomolecular systems ⁴. The simplest element is a gene whose protein $p$ is made at a rate proportional to an input $u(t)$ and removed by degradation and dilution at rate $\gamma$: $\frac{dp}{dt} = k\,u(t) - \gamma p$. Its transfer function is $\frac{P(s)}{U(s)} = \frac{k}{s + \gamma}$. To read this as a frequency response, substitute $s = i\omega$, where $\omega$ is the angular frequency of a sinusoidal input and $i$ is the imaginary unit. The magnitude of the response is $\frac{k}{\sqrt{\omega^2 + \gamma^2}}$. Inputs that vary slowly compared with the protein lifetime $1/\gamma$ pass through with gain $k/\gamma$. Inputs that vary faster are averaged away. A gene is a first-order low-pass filter whose bandwidth is set by how fast its protein disappears.

Circuits follow by algebra. A cascade of genes has a transfer function that is the product of the stages, so each stage adds a further low-pass filter and delay. A negative feedback loop with forward path $G(s)$ and feedback path $K(s)$ has closed-loop transfer function $\frac{G}{1 + GK}$. If the loop gain $GK$ is large the response approaches $1/K$, so the output is set by the feedback element and becomes insensitive to variation in the forward path. Del Vecchio, Dy and Qian review how these ideas have been carried into living cells and the practical obstacles met on the way, writing for both synthetic biologists and control engineers ⁵.

Measuring the frequency response

Rather than compute $H(s)$ from a model, one can measure it. Drive a cell with a periodic input, sweep the period, and record how much of the oscillation reaches the output. The magnitude and phase of the response as a function of frequency are the transfer function evaluated along $s = i\omega$. Three groups did this to yeast in 2008.

Mettetal and colleagues applied periodic osmotic stimuli to yeast and measured the frequency dependence of the osmo-adaptation pathway. From those measurements they identified a predictive model and found that the dynamics are dominated by fast negative feedback acting through the kinase Hog1 that does not require new protein synthesis. For large shocks a slower feedback through gene expression appears as well ⁶. Hersen and colleagues built a microfluidic device that switches the osmolyte concentration at chosen frequencies. The HOG pathway behaved as a low-pass filter: it integrates fast changes and tracks slow ones. The frequency data yielded kinetic rates that had not been measurable before and showed that the Ssk1 branch of the pathway integrates fast signals while the Ste11 branch does not ⁷. Bennett and colleagues did the same for a metabolic network by oscillating the carbon source between glucose and galactose. The galactose system was again a low-pass filter, and the response was faster than models predicted. The discrepancy was traced to galactose transcripts whose half-lives depend on the carbon source ⁸.

The caveat in all three is the one the transform always carries. A frequency response is a property of a linear system, and a cell is only linear for small enough perturbations around an operating point. The large-shock behaviour that Mettetal’s group observed is the pathway leaving that regime.

Generating functions: the discrete cousin

Mattuck’s power series returns in stochastic gene expression. When molecule counts are small, the state of a cell is a probability $P(n, t)$ of holding $n$ molecules at time $t$. The chemical master equation is the set of coupled differential equations that moves probability between neighbouring values of $n$ as molecules are made and destroyed. If $a(n) = P(n)$ is used as the coefficient sequence, the power series $G(x) = \sum_{n=0}^{\infty} P(n)\,x^n$ is the probability generating function. The master equation gives one equation per $n$, and together they become a single partial differential equation for $G(x, t)$. The moments come from derivatives of $G$ at $x = 1$, and if $G$ can be solved in closed form the whole distribution follows.

Thattai and van Oudenaarden analysed intrinsic noise in gene networks this way. They found that protein noise is controlled mainly at the level of translation, so that the mean and the variance of a protein can be tuned independently. Negative autoregulation reduces noise, and fluctuations in a regulatory protein propagate down a cascade ⁹. Shahrezaei and Swain took the two-stage model of an mRNA that is transcribed and translated and used the fact that in budding yeast proteins decay much more slowly than their mRNA. With that separation of timescales they derived the protein distribution in closed form, showed that protein is made in geometrically distributed bursts, and found the distributions to be asymmetric enough that mean and variance alone describe them poorly ¹⁰. The Laplace transform is the continuous limit of the same trick.

Survival analysis and frailty

The last example is statistics rather than dynamics. The survival function $S(t)$ is the probability that an event, such as death or relapse, has not occurred by time $t$. The hazard $h(t)$ is the instantaneous event rate among those still at risk, and its integral $H(t)$ is the cumulative hazard, with $S(t) = \exp(-H(t))$. A frailty model says that individuals differ by an unobserved positive multiplier $Z$ on the hazard. Someone with $Z = 2$ faces twice the baseline hazard at every age. Conditional on $Z$, survival is $\exp(-Z\,H(t))$.

The population survival curve is the average of that over the distribution of $Z$. Writing the expectation out, $S_{\text{pop}}(t) = E[\exp(-Z\,H(t))]$, and this is exactly the Laplace transform of the frailty distribution evaluated at the point $H(t)$. Hougaard’s monograph on multivariate survival data builds its treatment of frailty on this identity ¹¹. Aalen introduced a broad family of frailty distributions that contains the gamma distribution and also allows a non-susceptible fraction of the population, and applied it to breast cancer survival from the Norwegian Cancer Registry ¹². For gamma frailty with mean 1 and variance $1/\delta$ the transform is $(1 + u/\delta)^{-\delta}$, so $S_{\text{pop}}(t) = \left(1 + \frac{H(t)}{\delta}\right)^{-\delta}$. The consequence is selection: the frail leave the population first, so the observed population hazard flattens or declines even when every individual’s hazard is rising. Any study of a biomarker’s relation to an outcome over decades has to reckon with this.

What this means for your own time series

A glucose trace after a meal and a drug concentration after a dose are both responses to an input. The linear-systems view says the response is the input convolved with an impulse response. That impulse response is characterised by a handful of time constants: how fast the curve rises, when it peaks, and how fast it decays. We think this is the right first model for a person’s own data, for three reasons. Superposition means two meals close together should produce the sum of two single-meal responses, so a departure from that sum is information. Time invariance means the same meal on different days should produce the same curve, so a change in the curve is a change in the system rather than in the input. And residuals from a linear prediction are where the nonlinear physiology shows up.

The view breaks in known places. Glucose regulation is nonlinear at its core: in the minimal model that Bergman and colleagues selected from seven candidates, insulin action multiplies the glucose level rather than adding to it, and it was a nonlinear model that best predicted the data ¹³. Drug elimination saturates at high concentrations ³. Parameters drift with time of day and with sleep and exercise, so the system is only time-invariant over windows of hours. Large inputs push the system out of the linear regime, which is the same lesson as the large osmotic shocks in yeast. And sampling limits what can be identified. A continuous glucose monitor reporting every five minutes can resolve nothing faster than a ten-minute period, and a blood panel drawn monthly carries almost no dynamic information: those markers are levels to be tracked and compared. Where a curve looks abnormal, its interpretation belongs with a clinician.

Questions people also ask

What is the Laplace transform in plain terms? It is the continuous version of summing a power series. The coefficient sequence $a(n)$ becomes a function $f(t)$, the sum becomes an integral, and the base $x$ is renamed $e^{-s}$.

Why does the Laplace transform turn differential equations into algebra? Because differentiating in time corresponds to multiplying by $s$ in the transform domain, up to an initial-value term. Convolution of an input with an impulse response becomes a product.

Where is the Laplace transform used in biology? In compartmental pharmacokinetics to get concentration-time curves from dose inputs ¹ and in the control theory of gene circuits as transfer functions ⁴. It also appears in the frequency-response analysis of signaling pathways ⁶ and in survival analysis, where the transform of the frailty distribution gives the population survival curve ¹¹.

When does a linear-systems model of a biological response fail? When the response saturates or the input is large enough to leave the linear regime, when parameters change over the observation window, and when sampling is too sparse to identify the dynamics.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Milo Gibaldi, Donald Perrier. Method of Laplace Transforms. In Pharmacokinetics, Second Edition. Marcel Dekker (CRC Press), 1982. https://doi.org/10.1201/b14095-13 ↩ ↩²

2. Leslie Z. Benet. General Treatment of Linear Mammillary Models with Elimination from any Compartment as Used in Pharmacokinetics. Journal of Pharmaceutical Sciences, 1972. https://doi.org/10.1002/jps.2600610408 ↩

3. Milo Gibaldi, Donald Perrier. Nonlinear Pharmacokinetics. In Pharmacokinetics, Second Edition. Marcel Dekker (CRC Press), 1982. https://doi.org/10.1201/b14095-7 ↩ ↩²

4. Domitilla Del Vecchio, Richard M. Murray. Biomolecular Feedback Systems. Princeton University Press, 2015. https://doi.org/10.1515/9781400850501 ↩ ↩²

5. Domitilla Del Vecchio, Aaron J. Dy, Yili Qian. Control theory meets synthetic biology. Journal of The Royal Society Interface, 2016. https://doi.org/10.1098/rsif.2016.0380 ↩

6. Jerome T. Mettetal, Dale Muzzey, Carlos Gómez-Uribe, Alexander van Oudenaarden. The Frequency Dependence of Osmo-Adaptation in Saccharomyces cerevisiae. Science, 2008. https://doi.org/10.1126/science.1151582 ↩ ↩²

7. Pascal Hersen, Megan N. McClean, L. Mahadevan, Sharad Ramanathan. Signal processing by the HOG MAP kinase pathway. Proceedings of the National Academy of Sciences, 2008. https://doi.org/10.1073/pnas.0710770105 ↩

8. Matthew R. Bennett, Wyming Lee Pang, Natalie A. Ostroff, et al. Metabolic gene regulation in a dynamically changing environment. Nature, 2008. https://doi.org/10.1038/nature07211 ↩

9. Mukund Thattai, Alexander van Oudenaarden. Intrinsic noise in gene regulatory networks. Proceedings of the National Academy of Sciences, 2001. https://doi.org/10.1073/pnas.151588598 ↩

10. Vahid Shahrezaei, Peter S. Swain. Analytical distributions for stochastic gene expression. Proceedings of the National Academy of Sciences, 2008. https://doi.org/10.1073/pnas.0803850105 ↩

11. Philip Hougaard. Analysis of Multivariate Survival Data. Springer, 2000. https://doi.org/10.1007/978-1-4612-1304-8 ↩ ↩²

12. Odd O. Aalen. Heterogeneity in survival analysis. Statistics in Medicine, 1988. https://doi.org/10.1002/sim.4780071105 ↩

13. Richard N. Bergman, Y. Z. Ider, Charles R. Bowden, Claudio Cobelli. Quantitative estimation of insulin sensitivity. American Journal of Physiology-Endocrinology and Metabolism, 1979. https://doi.org/10.1152/ajpendo.1979.236.6.e667 ↩