Introduction to Mathematical Modelling¶

by Dr Kamil Erguler
Climate and One Health Group
CARE-C - The Cyprus Institute

---

Previous (Parameters) --- Index --- Next (Selection)

Thermal performance curves — Part 1¶

Temperature is an important driver of metabolism, especially in organisms whose body temperature depends strongly on the surrounding environment. Many insects, reptiles, amphibians, fish, and microorganisms fall into this broad category. They are often described as ectotherms, because they rely largely on external sources of heat, or poikilotherms, because their internal temperature can vary substantially.

One way to study this relationship is to measure a physiological process, such as growth, development, reproduction, or survival, at different constant temperatures. Each measurement gives us one point: a temperature and the corresponding level of performance.

Once we have such measurements, we can try to describe them with a mathematical equation. This allows us to summarise the observations, fill in values between measured temperatures, and make predictions under new conditions.

Corkrey and colleagues published a comprehensive analysis of temperature–growth relationships in 2016. Let's examine a small part of their data for the common gut bacterium, Escherichia coli, and see what kind of model we can build from it.

From this plot, it seems that there may be a roughly linear relationship between bacterial growth performance and temperature. When it is cold, E. coli grows slowly, and therefore food may take longer to spoil. This is comforting in winter, but you should not leave food on the counter for long in summer.

Let us propose a model for this relationship, and call it $M$.

Based on what we have seen so far, we may describe the temperature–performance relationship by

$$ \rho(T) = aT + b, $$

where $\rho(T)$ is the temperature-dependent growth rate, $T$ is temperature, $a$ controls how strongly the growth rate increases with temperature, and $b$ is an offset.

Let us restrict $a$ to non-negative values. This means that $M$ is not just any linear model; it is the positive-association version of a linear relationship model structure. A separate model version could instead focus on the negative domain of $a$, representing the hypothesis that growth rate decreases as temperature increases. Here, to avoid complicating things too early, we will consider only the positive-association case.

---

An obvious limitation is that this model describes a straight line: one predicted performance value for each temperature, with performance changing at a constant rate as temperature changes. From the plot, it is clear that no single straight line can pass through all the points.

Why not? Experimental error? Differences in conditions other than temperature? Natural variability among replicate measurements?

It could be many things.

So, we need to account for variability in the observations. Otherwise, the model cannot match the data. We can introduce this variability with a noise term from a chosen probability distribution:

$$ \rho(T) = aT + b + \epsilon, $$

where $\epsilon$ is a random number. Let us say that it is normally distributed around zero, with variance $\sigma^2$:

$$ \epsilon \sim N(0, \sigma^2). $$

Therefore, $M$ has three parameters, $a$, $b$, and $\sigma$, and gives a distribution of possible growth rates around the straight-line relationship.

It seems that, with $a = 0.004$, $b = -0.08$, and $\sigma = 0.01$, our model does not quite match the observations. We can quantify this mismatch using the likelihood of the observations, $\delta$, defined as the probability density assigned to the data by the model under a given parameter configuration.

Strictly speaking, for discrete observations, the likelihood is the probability assigned to the data. For continuous observations, the likelihood is the probability density assigned to the data. Since growth rate is treated as a continuous variable in this example, we calculate the likelihood using a probability density:

$$ P(\delta \mid M, \theta). $$

For any given temperature, $T$, the model generates a random growth rate with mean $aT + b$ and standard deviation $\sigma$. Now consider one data point, where the temperature is $T_i$ and the measured growth rate is $\delta_i$. The likelihood of this observation is the probability density that the model assigns to $\delta_i$ at temperature $T_i$:

$$ P(\delta_i \mid T_i, M, \theta) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[ -\frac{(\delta_i - (aT_i+b))^2}{2\sigma^2} \right]. $$

Assuming that the observations are independent of each other, the likelihood of the full dataset is the product of these probability densities:

$$ P(\delta \mid M, \theta) = \prod_i \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[ -\frac{(\delta_i - (aT_i+b))^2}{2\sigma^2} \right]. $$

A note on probability density¶

You may notice that the likelihood can sometimes be a very large number. This does not mean that the model assigns a very large probability to the observations. In this case, the likelihood is calculated using probability densities, and probability densities are not probabilities. For continuous random variables, probabilities are assigned to intervals, not to exact values.

To see why, imagine drawing a random number from a uniform distribution between 0 and 1. The probability density is the same at 0.5 and at 0.5001. The probability of drawing a number within the interval $[0.49999, 0.50001]$ is also the same as the probability of drawing a number within $[0.50009, 0.50011]$, because the two intervals have the same width. However, the probability of drawing exactly 0.5 is zero. To include only 0.5 and no other real number, the interval would have to have zero width.

Setting up priors¶

How do we quantify our prior belief about the three parameters?

$$ \theta = (a,b,\sigma) $$

Is there a reason to lean towards a certain region of parameter space?

Testing a specific hypothesis can be a reason to do so. For instance, we restricted $a$ to non-negative values to represent the hypothesis that growth rate does not decrease as temperature increases.

We must also restrict $\sigma$ to positive values, because it represents the standard deviation of the noise term, and a standard deviation must be greater than zero.

The easiest way to express these restrictions is to assign uniform prior distributions with predefined ranges. Let us choose the following:

$$ \begin{array}{lcl} P(a) &=& U(a \mid \min=0,\max=0.01),\\ P(b) &=& U(b \mid \min=-0.2,\max=0.2),\\ P(\sigma) &=& U(\sigma \mid \min=0.0001,\max=0.10). \end{array} $$

These priors do not favour any particular value within the allowed ranges. They simply define the part of parameter space that this model is allowed to explore.

Using the likelihood and the prior, we can define the unnormalised posterior density:

$$ P(\theta \mid \delta, M) \propto P(\delta \mid \theta, M)P(\theta \mid M). $$

In the code, we will work with the log of the unnormalised posterior for numerical stability:

The posterior density plot shows which combinations of $a$ and $b$ are most supported by the observations, while keeping $\sigma = 0.01$ fixed. Each point in the $(a,b)$ parameter grid defines one possible model relationship:

$$ \rho(T) = aT + b + \epsilon. $$

However, not all parameter combinations are equally supported by the observations. Within the allowed parameter range, we observe a strong dependence between $a$ and $b$. Larger values of $a$ must be paired with lower values of $b$ to keep the line close to the observed data. In other words, $a$ and $b$ show a negative association in the posterior.

The negative values of $b$ should be interpreted cautiously. The intercept is the predicted growth rate at $T=0$, which is outside the temperature range shown here. In this example, $b$ mainly helps position the line within the observed temperature range, rather than carrying a direct biological meaning.

Markov chain Monte Carlo (MCMC)¶

So far, we have been quite lucky. We were able to search through the interesting parts of parameter space, identify the posterior mode, and calculate the expected relationship between temperature and growth rate, while still keeping $\sigma$ fixed.

MCMC offers a powerful alternative. Instead of evaluating the posterior over a predefined grid, we allow the algorithm to move through parameter space step by step. After an initial exploration period, the algorithm generates a sequence of parameter configurations that represents the posterior distribution. Regions with higher posterior density are visited more often, while lower-density regions are visited less often.

The Metropolis–Hastings algorithm¶

Metropolis–Hastings is an MCMC algorithm based on proposals and acceptances. At each step, it proposes a new parameter configuration near the current one, then decides whether to move there by comparing the current and proposed posterior support.

This accept-or-reject rule turns an ordinary random walk into a sampler for the posterior distribution.

Here, let's begin the chain at an initial parameter configuration, for example

$$ \theta = (a,b,\sigma) = (0.004, -0.08, 0.01), $$

At each step, the algorithm proposes a small random change to the current parameter values,

$$ \theta' = (a',b',\sigma'). $$

This proposal is drawn from a proposal distribution, which we can write as

$$ q(\theta' \mid \theta). $$

This means: given the current configuration $\theta$, the proposal distribution tells us how likely we are to propose the new configuration $\theta'$.

We then compare the posterior value of the proposed configuration with the posterior value of the current configuration. If the proposal is accepted, we move to the new position and make another proposal from there. If it is rejected, we stay where we are and try again.

---

Whether the proposal is accepted or not depends on the relative posterior values of the current configuration, $P(\theta \mid \delta, M)$, and the proposed configuration, $P(\theta' \mid \delta, M)$, as well as the relative probabilities of proposing $\theta'$ from $\theta$ and proposing $\theta$ from $\theta'$.

The posterior part is intuitive: we want the sampler to spend more time in high-density regions of the posterior. Therefore, proposals that increase the posterior value are more likely to be accepted.

The proposal part is a little trickier. If it is easier to propose the move $\theta \rightarrow \theta'$ than the reverse move $\theta' \rightarrow \theta$, then accepting the move too readily would bias the sampler: it would spend too much or too little time in some parts of parameter space simply because of how proposals are generated. The Metropolis–Hastings correction offsets this imbalance.

---

We can formulate the above criteria into an acceptance ratio, written as

$$ R = \frac{ P(\theta' \mid \delta, M) }{ P(\theta \mid \delta, M) } \frac{ q(\theta \mid \theta') }{ q(\theta' \mid \theta) }. $$

This is equivalent to defining an acceptance probability as

$$ \alpha = \min(1, R). $$

When $R > 1$, the proposed configuration is always accepted. When $R < 1$, the proposed configuration is accepted with probability $R$. For example, if the posterior value at the proposed configuration is half as large as the posterior value at the current configuration, and the proposal distribution is symmetric, then $R = 0.5$, $\alpha = 0.5$, and the proposal is accepted with probability 50%.

This is the part that prevents the sampler from only climbing uphill. Occasionally accepting worse parameter configurations allows it to explore the surrounding parameter space instead of getting stuck near the first good configuration it finds.

After many steps, the sampler spends more time in regions of parameter space with higher posterior probability. The resulting sequence of parameter values is an approximation of the posterior distribution.

Working with ratios¶

Using Bayes' rule, we can expand the posterior terms in the acceptance ratio:

$$ R = \frac{ P(\theta' \mid \delta, M) }{ P(\theta \mid \delta, M) } \frac{ q(\theta \mid \theta') }{ q(\theta' \mid \theta) }. $$

For the current configuration,

$$ P(\theta \mid \delta, M) = \frac{ P(\delta \mid \theta, M)P(\theta \mid M) }{ P(\delta \mid M) }, $$

and for the proposed configuration,

$$ P(\theta' \mid \delta, M) = \frac{ P(\delta \mid \theta', M)P(\theta' \mid M) }{ P(\delta \mid M) }. $$

Substituting these into the acceptance ratio gives

$$ R = \frac{ \frac{ P(\delta \mid \theta', M)P(\theta' \mid M) }{ P(\delta \mid M) } }{ \frac{ P(\delta \mid \theta, M)P(\theta \mid M) }{ P(\delta \mid M) } } \frac{ q(\theta \mid \theta') }{ q(\theta' \mid \theta) }. $$

The model evidence, $P(\delta \mid M)$, appears in both the numerator and the denominator, so it cancels out:

$$ R = \frac{ P(\delta \mid \theta', M)P(\theta' \mid M) }{ P(\delta \mid \theta, M)P(\theta \mid M) } \frac{ q(\theta \mid \theta') }{ q(\theta' \mid \theta) }. $$

This is useful because we do not need to calculate the evidence explicitly. We only need the likelihood, the prior, and the proposal probabilities.

If the proposal distribution is symmetric, such as a normal random-walk proposal with fixed standard deviations, then the move $\theta \rightarrow \theta'$ and the reverse move $\theta' \rightarrow \theta$ are equally likely:

$$ q(\theta \mid \theta') = q(\theta' \mid \theta). $$

In that case, the proposal terms cancel, and the acceptance ratio can be written as

$$ R = \frac{ P(\delta \mid \theta', M)P(\theta' \mid M) }{ P(\delta \mid \theta, M)P(\theta \mid M) }. $$

This is the ratio of the unnormalised posterior densities we calculated earlier.

If the prior is uniform, and both the current and proposed configurations are inside the allowed parameter ranges, then

$$ P(\theta' \mid M) = P(\theta \mid M), $$

so the prior terms also cancel. In that case, we are left with the ratio of the likelihoods:

$$ R = \frac{ P(\delta \mid \theta', M) }{ P(\delta \mid \theta, M) }. $$

Therefore, with a symmetric proposal distribution and uniform priors, Metropolis–Hastings accepts or rejects proposed moves by comparing how much likelihood the proposed and current parameter configurations assign to the observations.

Burn-in and thinning¶

The Metropolis–Hastings sampler generates a sequence of parameter configurations, called a Markov chain. However, we should not necessarily treat every sampled value as an independent draw from the posterior distribution.

Two practical issues are important.

First, the chain starts from an initial parameter configuration. Early in the run, the sampler may still be moving away from this starting point and towards the high-density region of the posterior. These early samples can therefore depend strongly on the initial values. To reduce this effect, we usually remove an initial part of the chain. This discarded part is called the burn-in. A common rule of thumb is to discard the first quarter of the output, although in practice this should depend on how quickly the chain appears to settle.

Second, consecutive samples in an MCMC chain are not independent. Each proposal is made from the current position, so nearby iterations tend to be correlated with each other. This is called autocorrelation. If autocorrelation is high, the chain may contain many samples, but less independent information than the number of rows suggests.

One simple way to reduce autocorrelation among retained samples is thinning. This means keeping only every $k$-th sample from the chain. For example, with $k=10$, we keep every tenth sample and discard the rest.

More formal MCMC analyses use diagnostics such as trace plots, autocorrelation plots, effective sample size, and comparisons among multiple chains started from different initial values.

Posterior predictive distribution¶

The posterior predictive distribution of our model, $M$, is obtained by averaging the model predictions over all possible parameter configurations, weighted by the posterior support for each configuration.

This can be understood as adding another layer of variability to an already stochastic model. Our model $M$ is stochastic because it describes a distribution of possible growth rates, rather than a single growth rate. It does this using three parameters: $a$ and $b$, which describe the linear relationship with temperature, and $\sigma$, which describes the variability around this relationship.

The model would have been simpler if it were deterministic: that is, if it always generated the same growth rate for the same temperature. But then it would not explain the variability we see in the data.

For each parameter configuration,

$$ \theta_j = (a_j, b_j, \sigma_j), $$

we obtain a different linear relationship and a different probability distribution for growth rate:

$$ \rho_j(T) = a_jT + b_j + \epsilon_j, $$

where

$$ \epsilon_j \sim N(0,\sigma_j^2). $$

The posterior predictive distribution is defined by the integral

$$ P(\rho(T) \mid T, \delta, M) = \int P(\rho(T) \mid T, \theta, M) P(\theta \mid \delta, M) \,d\theta . $$

This integral averages the model predictions over all possible parameter configurations, weighted by their posterior support, $P(\theta \mid \delta, M)$. As a result, we obtain a single predictive distribution supported by the data, rather than a prediction based on one particular parameter configuration.

Monte Carlo integration¶

We can approximate this integral using samples from the MCMC chain.

If we have $J$ posterior samples,

$$ \theta_j = (a_j,b_j,\sigma_j), $$

then the posterior predictive density can be approximated by

$$ P(\rho(T) \mid T, \delta, M) \approx \frac{1}{J} \sum_{j=1}^{J} P(\rho(T) \mid T, \theta_j, M). $$

For our model,

$$ P(\rho(T) \mid T, \theta_j, M) = N(\rho(T) \mid a_jT+b_j, \sigma_j^2). $$

So, at each temperature, we calculate the normal density implied by every posterior sample and then average these densities. This is Monte Carlo integration: the integral over parameter space is replaced by an average over sampled parameter values.

To illustrate this, we plot the posterior predictive density at $T=25^{\circ}\mathrm{C}$. This density is the average of many normal densities, one from each posterior sample. In other words, it is a mixture of normals.

Is the model predictive?¶

Within the measured temperature range, the model performs reasonably well. This means that we can estimate, with reasonable accuracy, what the growth rate might be for any temperature between $15^{\circ}C$ and $35^{\circ}C$. Within this range, the model assumes that growth rate changes linearly with temperature.

But what about temperatures outside this range? Can we also predict growth rate at $10^{\circ}C$ with this model? It appears that we cannot. As we see from the figure, the expected growth rate falls below zero at temperatures lower than $15^{\circ}C$ according to the model. An eerie prediction indeed.

We could say that we are not interested in temperatures outside the observed range and use the linear model only within the range of measured temperatures. This is called interpolation.

Using the model outside the observed range is called extrapolation. Extrapolation means asking the model to make predictions in regions where we have no observations. This can be useful for exploration, but it should be treated with great caution.

In this case, restricting the linear model to interpolation is a practical admission that the model is not meant to explain the full biological system. It may still be useful within the observed range, but we should not assume that the same straight-line relationship continues beyond it.

So, what should we do? The best advice is to collect more data and then return to improving the model, and our understanding of the system, based on this newfound knowledge.