The Wikipedia article on Diatoms states that:
an assemblage of living diatoms doubles approximately every 24 hours by asexual multiple fission; the maximum life span of individual cells is about six days
So a simple exponential equation doesn't suffice anymore to describe population growth over time, since dead organisms don't contribute to it anymore.
Is there a formula that can take this limited lifespan into account (ignoring all other possible constraints like space or energy)?
As a first approximation, we can use a linear discrete time model with $t \in \mathbb{N}$ being days:
Let's assume that our population can be partitioned into seven states based on age (in days) $S:=\{0,1,2,3,4,5,D\}$ where $0$ is the number of diatoms that were just "born" and $D$ is just a placeholder for all diatoms 6 days and older (i.e., dead).
Note I am not accounting for any distribution in lifespans (e.g., using an exponential distribution of lifetimes) -- which assumes that the range is rather constrained (e.g., rare for a diatom to die at day 3 or live past day 7).
Let $P_t \in \mathbb{N^7}$ be the vector of population counts by age at day $t$. Let's assume that when we go from $t$ to $t+1$ that diatoms undergo binary fission, with one half of the fission generating "newborns" and the other half leaving "aged" diatoms (one day older).
This means that if I have just $2$ newborn diatoms at time $0$ (so $P_0=(2,0,0,0,0,0,0)$) then I will have $P_1=(2,2,0,0,0,0,0)$ at $t=1$ and $P_2=(4,2,2,0,0,0,0)$ at $t=1$ etc. with the exception that any diatoms in state $D$ (i.e., $(0,0,0,0,0,0,N)$) stay there and don't add births.
We can represent this process as a linear transformation where our initial population vector $P_0$ is left-multiplied by a population dynamics matrix $A$:
$$A = \begin{bmatrix}1&1&1&1&1&1&0\\1&0&0&0&0&0&0\\0&1&0&0&0&0&0\\0&0&1&0&0&0&0\\0&0&0&1&0&0&0\\0&0&0&0&1&0&0\\0&0&0&0&0&1&1 \end{bmatrix}$$
With the above definition for $A$, we can get $P_t$ as follows:
$$P_t = A^tP_0$$
Unfortunately, this is only partly helpful, because the population will grow without bound. But we can also get the population distribution by normalizing by the total population at time $t$ (i.e., the $L_1$ norm of $P_t$):
$$p_t = \frac{P_t}{||P_t||_1}$$
There may be some nice formula here but I went the easy route and coded this up using Python + sympy (you could use numpy too but sympy is easier IMO for some stuff):
from sympy.matrices import Matrix
top_row = [1]6+[0]
mid_rows = [[0]i + [1] + [0](6-i) for i in range(5)]
last_row = [0]5 +[1,1]
rows = [top_row]+ mid_rows + [last_row]
A = Matrix(rows)
def popcounts(A,n,y_0=None):
if not y_0:
y_0=Matrix([1]7)
return (An)y_0
def popdist(A,n,y_0=None):
Pn = popcounts(A,n,y_0)
return (1/float(sum(Pn)))*Pn
What you'll find is that popdist(A,600) is invariant to the choice of $P_0$ and converges (numerically) to the following (rounding to 2 decimal places):
$P_{\infty} = (0.50, 0.25, 0.13, 0.05, 0.03, 0.02, 0.02)$
So with this simple linear model, we get that the population converges to 50% newborns and the majority being 0, 1 or 2 days old.
The fact that there are 25x more newborns than total dead shows that the population is unbounded and that this model clearly has limitations (e.g., they live 6 days assuming ample nutrients but lifespan drops rapidly as deprivation takes hold).
I'm sure we can think of many refinements to this model (and I'm sure biologists have better models already) but this is the first approximation I could think of.
So I gave this a little thought and realized that we can just use the eigen-decomposition of $A$ to get more insight into the limiting behavior of $A^n$
The eigendecomposition of A involves the matrix of eigenvectors $Q$ and the associated diagonal matrix $D$ where the entries are the eigenvalues of the associated eigenvectors in $Q$ (same order in both).
We can represent $A$ using $Q,D$ as follows:
$$A = QDQ^{-1} \implies A^n = QD^nQ^{-1}$$
So matrix powers are ridiculously easy if the eigendecomposition exists (doesn't always exist -- $A$ must be diagonalizable)
What we'll do is look at the largest absolute eigenvalue $\lambda^* = \max_i |\lambda_i|$ to see what will happen as $n\to \infty$.
In general, we can have:
**Case 1: ** $\lambda^*<1 \implies ||A^n||_{\infty}\to \mathbf{0}$
**Case 2: **$\lambda^*=1 \implies ||A^n||_{\infty} \to c < \infty $
**Case 3: ** $\lambda^*> 1 \implies ||A^n||_{\infty} \to \infty$
In Case 2 we have $D^{n}$ converging to a diagonal matrix $D_{\infty}$ where all eigenvalues $< 1$ are set to zero and the eigenvalue(s) equal to $1$ are left as is.
Therefore,
$$A_{\infty} = QD_{\infty}Q^{-1}$$
In Case 3 $D_{\infty}$ will not exist since at least one entry will grow without bound and dominate the others.
However, perhaps we can argue that a normalized version of $A_{\infty}$ (let's call it $\alpha$) does exist if we replace $D$ with a normed version:
$$\Delta_n:=\frac{1}{\text{tr}\left(D^n\right)}D^{n}$$
Then we have:
$$\text{tr}\left(D^n\right) \sim O((\lambda^*)^n) \implies \Delta_n \to Z(i)$$
Where $Z(i)$ is a matrix that has zeroes everywhere except $Z_{i,i} = 1$, where $i = \arg \max_i |\lambda_i|$.
Therefore,
$$\alpha = QZ(i)Q^{-1}$$
If our starting population state $P_0$ is as before (uniform across all states), then the limiting "normed" population vector $P_{\infty}$ is:
$$P_{\infty} = \alpha p_0$$
However, this not a distribution, but a limiting of vector of states that reflects the asymptotic relative sizes of the states as $n\to \infty$.
To get our long-run distribution $p_{\infty}$ we can divide $P_{\infty}$ by it's $L_1$ norm to get it so sum to 1:
$$p_{\infty} = \frac{P_{\infty}}{||P_{\infty}||_1}$$
I wrote up a small python script to check this logic (not a proof obviously):
>>> import numpy.linalg as la
>>> import numpy as np
>>>
>>> # set up our model matrix
>>> top_row = [1]*6+[0]
>>> mid_rows = [[0]*i + [1] + [0]*(6-i) for i in range(5)]
>>> last_row = [0]*5 +[1,1]
>>> rows = [top_row]+mid_rows + [last_row]
>>> A = np.array(rows)
>>> print("A is \n", A,"\n", "with dimensions: ", A.shape)
A is
[[1 1 1 1 1 1 0]
[1 0 0 0 0 0 0]
[0 1 0 0 0 0 0]
[0 0 1 0 0 0 0]
[0 0 0 1 0 0 0]
[0 0 0 0 1 0 0]
[0 0 0 0 0 1 1]]
with dimensions: (7, 7)
>>>
>>> # Let's look at the size of the eigenvalues
>>> D,U=la.eig(A)
>>> print("The absolute value of the eigenvalues are ",np.round(np.absolute(D),2))
The absolute value of the eigenvalues are [1. 1.98 0.84 0.85 0.85 0.91 0.91]
>>> # we see that the second eigenvalue is the only one greater than 1
>>> # therefore, the second eigenvalue will dominate as we take greater powers of A
>>>
>>> # We're going to calculate what the (normalized) limiting matrix will be
>>> # We'll do this by noticing that (1/tr(|D^n|))D^n} approaches a matrix with all zeros except at D[1,1] where it equals 1
>>> # Let's calculate this limiting "normed" D
>>> d = np.zeros((7,7))
>>> d[1,1]=1
>>> print("Normed diagonal limit is \n", d)
Normed diagonal limit is
[[0. 0. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]]
>>>
>>> # Now we can get the limiting (normed) verison of A
>>> A_lim = U@d@la.inv(U)
>>> print("A_lim is \n", np.round(A_lim,2))
A_lim is
[[0.52+0.j 0.51+0.j 0.5 -0.j 0.46+0.j 0.4 +0.j 0.26-0.j 0. +0.j]
[0.26+0.j 0.26+0.j 0.25-0.j 0.23+0.j 0.2 +0.j 0.13-0.j 0. +0.j]
[0.13+0.j 0.13+0.j 0.13-0.j 0.12+0.j 0.1 +0.j 0.07-0.j 0. +0.j]
[0.07+0.j 0.07+0.j 0.06-0.j 0.06+0.j 0.05+0.j 0.03-0.j 0. +0.j]
[0.03+0.j 0.03+0.j 0.03-0.j 0.03+0.j 0.03+0.j 0.02-0.j 0. +0.j]
[0.02+0.j 0.02+0.j 0.02-0.j 0.02+0.j 0.01+0.j 0.01-0.j 0. +0.j]
[0.02+0.j 0.02+0.j 0.02-0.j 0.02+0.j 0.01+0.j 0.01-0.j 0. +0.j]]
>>>
>>> # Now we'll set up our starting distribution (where we assumed uniform aross all states)
>>> p0 = np.ones((7,1))
>>>
>>> # Finally, we'll calcualte the limting (normed) population
>>> p_lim = A_lim@p0
>>> print("Normed limiting population is:\n", np.round(np.absolute(p_lim),3))
Normed limiting population is:
[[2.652]
[1.337]
[0.674]
[0.34 ]
[0.171]
[0.086]
[0.088]]
>>>
>>> # The final step is to convert the normed population to a distribution
>>> dist_lim = np.absolute(p_lim*(1/la.norm(p_lim,ord=1)))
>>> print("The limiting distribution of states is:\n", np.round(dist_lim,2))
The limiting distribution of states is:
[[0.5 ]
[0.25]
[0.13]
[0.06]
[0.03]
[0.02]
[0.02]]
>>>
Which agrees to within rounding error.