On models of physiologically structured populations and their reduction to ordinary differential equations.

Considering the environmental condition as a given function of time, we formulate a physiologically structured population model as a linear non-autonomous integral equation for the, in general distributed, population level birth rate. We take this renewal equation as the starting point for addressing the following question: When does a physiologically structured population model allow reduction to an ODE without loss of relevant information? We formulate a precise condition for models in which the state of individuals changes deterministically, that is, according to an ODE. Specialising to a one-dimensional individual state, like size, we present various sufficient conditions in terms of individual growth-, death-, and reproduction rates, giving special attention to cell fission into two equal parts and to the catalogue derived in an other paper of ours (submitted). We also show how to derive an ODE system describing the asymptotic large time behaviour of the population when growth, death and reproduction all depend on the environmental condition through a common factor (so for a very strict form of physiological age).

Introduction

Monod (1942, 1949, 1950; Morin and Monod 1942) presented his now famous ODE-model of batch and continuous cultivation of bacteria feeding on a (growth limiting) substrate:where X(t) is the bacterial biomass and S(t) is the substrate concentration in the growth tank at time t. The constant concentration of substrate in an external reservoir is denoted by and this is fed into the growth tank at a constant rate D. Bacteria and substrate are removed at the same rate, which therefore is called the dilution rate. The reciprocal of is a conversion factor telling how consumed substrate is converted into bacterial biomass. The case corresponds to batch cultivation.

The model is formulated in terms of bacterial biomass and not in terms of individual cells. The reason is that consumption of substrate results in cell growth and not directly in cell fission. When a cell splits depends on its size, which in turn depends on the history of its substrate consumption, so on the history of the substrate concentration and uptake ability as determined by size. The model should therefore be viewed as a reduction of a more detailed structured population model in the spirit of Tsuchiya et al. (1966), Fredrickson et al. (1967), Gyllenberg (1982) and Diekmann et al. (1984, 1986). The following question arises naturally: Under what conditions on the individual processes, viz. growth, survival and fission, do the Eqs. (1.1) and (1.2) provide a faithful representation of a size-structured cell-based model?

In two recent papers (Diekmann et al. 2018, submitted) we considered the problem of ODE-reducibility of physiologically structured populations from different angles and under different restrictions. In the paper (Diekmann et al. 2018) we focussed on the population level birth rate under the assumption of a finite number of states-at-birth and with individual state (i-state for short) dynamics, for instance, consisting of jumps from one state to another. We gave necessary and sufficient conditions on the kernel of the resulting renewal equation with input for when the birth rate vector can be recovered from the solution of a system of finitely many ODEs. In the paper (Diekmann et al. submitted) we considered deterministic i-state development and, working in the spirit of abstract evolution equations, we gave a complete catalogue, in terms of the individual growth and death rates of models, for which one can find a finite set of population outputs that can be recovered from a system of ODEs (see Sect. 6). The restrictions imposed were: the individual state space is one-dimensional and reproduction is part of the output.

The aim of the present paper is

to promote the renewal equation formulation of physiologically structured population models with deterministic i-state development;

to use that formulation to give a short and easy derivation of the key conditions for ODE reducibility;

to provide new examples, in particular for reproduction by fission into two equal parts;

to recall an old example (Diekmann et al. 1983; Metz and Diekmann 1986, I.4.5, II.14) of asymptotic ODE reducibility and to show how this fits in with the key condition.

Physiologically structured population models formulated as delay equations

Traditionally physiologically structured population models are formulated as first order partial differential equations with non-local terms (Metz and Diekmann 1986; Perthame 2007). In the papers (Diekmann et al. 1998, 2001) we developed a constructive method for solving such equations via generation expansion as it is used in the theory of renewal equations (Feller 1971). In the paper (Diekmann et al. 2010) we put the renewal equation centre stage by using the history of the birth rate as state and by defining a dynamical system by updating this history. In the present paper we take the renewal equation as the starting point for modelling physiologically structured populations.

Let be the individual state space (i-state space for short), that is, the set of all admissible individual states (i-states). The basic ingredient of a structured population model iswhich is the rate at which offspring are produced at time t in the Borel measurable set by an individual who had i-state at time s, given the environmental input E as a function of time. Note that Diekmann et al. (1998, 2001) considered the more general class of models that allow for jumps in the cumulative number of offspring produced as a function of the time elapsed, and therefore worked with the Stieltjes integral with respect to the time variable.

More often than not the possible states-at-birth form a subset which is considerably smaller than . As an extreme, but important, special case we have the situation of only one state-at-birth when . We therefore consider as a Borel measure on with support in . One of the reasons for working with measures instead of densities is that we can treat the cases of continuously and discretely distributed states-at-birth, as well as more complicated cases, in one go. We denote the -algebra of Borel subsets of a subset Y of by .

Let b(t) be the population birth rate, that is, b(t) is a Borel measure on with the interpretation that is the rate at which offspring with i-state in the set are produced at time t. Straightforward book-keeping then yields the following renewal equation with input:(Diekmann et al. 2007). If the birth rate is given on the interval as an initial conditionthen (2.1) and (2.2) can be written aswhereis the contribution to the birth rate of the individuals born before the initial time s.

Equation (2.1) can be viewed as a linear non-autonomous Volterra integral equation for functions of time with values in the Banach space of all Borel measures on . Results on existence and uniqueness of such abstract equations are available. The proofs usually mimic the corresponding ones for equations in (Gripenberg et al. 1990) based on Banach’s fixed point theorem and the associated Picard–Lindelöf successive approximations, see e.g. Kvapiš (1967). In our case well-posedness is almost immediate, because the solution can be constructed by way of the generation expansion (Diekmann et al. 1998, 2001).

We now introduce a second ingredient which is the probability that an individual that had i-state at time s is still alive at time t and has i-state in at that time, given the course of the environment (Diekmann et al. 1998, 2001). Once the birth rate b(t) has been solved from (2.1) the i-state distribution , often called the population state or p-state, can be obtained from the now explicit formulaWe now specialise to the case of deterministic i-state development which we simply call growth. Let g(x, E) be the individual growth rate of an individual with i-state x when the environmental condition is E. Then the i-state at time t of an individual, who had i-state at time s, equals x(t), where x is the solution to the initial value problemWe denote the rate at which individuals of i-state x produce children with i-state at birth in the measurable set , when the environmental condition is E, by and the death rate by . The probability that an individual, who had i-state at time s, is still alive at time t, given the course of the environmental condition (input), isIn terms of and we haveand hence for deterministic i-state development (2.1) and (2.5) can be written asIf the environmental input is a given function of time, then all information is contained in the linear non-autonomous integral equation (2.10) for the birth rate b. But individuals affect their own environment, for instance by consuming food and exploiting other resources, in other words, there is feedback to the environment. This feedback mechanism is modelled by integro-differential equations prescribing how the histories of b and E affect the current rate of change of E (Diekmann et al. 2010).

ODE reducibility

By a population output we mean a linear operator depending on the prevailing value E of the environmental input. Here we concentrate on outputs of the formwherewith w a measurable function from to .

In many concrete models the output weight function depends on the environment and so one should replace w in (3.2) by . However, it turns out that then an ODE representation is possible only ifwhich brings us back to (3.1) and (3.2).

If the population dynamics follows Eqs. (2.10) and (2.11), then the time evolution of the vector N(t) is given byWe want to derive a condition for when the population output N(t) satisfies an ordinary differential equationfor some -matrix valued function K of the environmental condition. Differentiating (3.4) and taking (2.10) into account, one finds that (3.5) holds ifwhere Dw(x) is the derivative of w at x represented by the (Jacobian) matrix. In the case of one-dimensional i-state space (), (3.6) of course simplifies toIn most of the rest of the paper we restrict to the case .

As explained in Sect. 2 one has, whenever feedback is present, to supplement the equation for the birth rate by an integro-differential equation describing the time evolution of the environment. Now, when we are looking for ODE reductions, this equation modelling the feedback should, of course, also be an ODE. If we know what properties of individuals are involved in the feedback, we know how to choose w and (3.7) becomes a condition on the model ingredients and . It may very well happen that (3.7) does not hold if we take as components of w the functions that we really need to describe the feedback, but that (3.7) does hold if we add additional components (Diekmann et al., submitted, Example 1.2). A systematic procedure for finding such additional components is to compute the left-hand side of (3.7), eliminate everything that is a linear combination (with E-dependent coefficients) of the current components, see whether what remains can be written as a linear combination (again with E-dependent coefficients) of some new functions of x and then add these new functions as additional components. Now repeat. This procedure is described by Diekmann et al. (submitted) in more detail as a generally applicable test for ODE reducibility starting from a collection of output functionals that one needs or wants.

Alternatively one might try to characterise combinations of and w for which one can find K(E) such that (3.7) holds. In the paper (Diekmann et al., submitted) we restricted the submodel for reproduction by requiring thatfor some -matrix M(E). Note that this does indeed hold if has the formwhich motivates us to describe this situation as “reproduction is part of the output”. We then found a complete catalogue of combinations of and w, parametrised by k and various functions of either x or E (see and in Sect. 6).

In the present paper we extend the material of (Diekmann et al., submitted) by providing some examples where (3.8) does not hold and by describing how large time ODE reduction is possible for models in which all physiological processes scale with the same factor depending on the environmental condition, so for a very strict form of physiological age.

One-dimensional output and Monod’s model

When the population output is one-dimensional, that is, when w is a scalar function (), the condition (3.7) says that the expressionshould be independent of x.

Next we apply this result to the Monod model discussed in the introduction. To conform with the notation used there, we now denote the environmental condition by S and interpret it as substrate. The Monod model is formulated in terms of biomass X (which now plays the role of N), soBacteria reproduce by fission into equal parts which we assume to happen at a rate depending on the size x of the mother cell and the substrate concentration S. We further assume that fission is successful with probability p(x, S) in which case the two daughter cells are both half the size of their mother; otherwise the mother cell is lost without leaving any offspring. We neglect any other kind of mortality. These assumptions mean that the reproduction rate is given byWith (4.2) and (4.3) we getAccording to our assumptions, cell loss happens at fission when the mother cell is lost at rate and through washout at rate D. The total loss rate is the sum of these rates:The condition for reducibility to the ODE (1.1) then becomesIn the case that every division is successful and (4.6) becomesthat is, individual growth rate per individual mass is exactly equal to the population growth rate per total biomass. The fission rate does not affect the finite dimensional representation in any way and can be chosen arbitrarily. When division may be unsuccessful, this type of mortality has to be compensated for in the individual growth rate and we haveNote that in this case the individual growth rate and the fission rate are related.

Two-dimensional output: a general formula and a particular example

We now describe how the previous class of examples can be extended. Our long term goal is to find a “maximal extension” in a similar manner as we did for the model class described in Sect. 6 below.

If , one can interpret (3.7) as two linear equations in the two unknowns g and . With the notational conventionand while suppressing both arguments x and E, these equations areWe now solve (5.2) in a somewhat peculiar way in order to tie up with in Sect. 6. Since(note that and need to be linearly independent functions of the variable x so the Wronskian vanishes at most in isolated points), we obtainwhere v is defined byIf reproduction is by division into two equal halves and if we incorporate the disappearance of the mother in rather than in , we haveNext we chooseso if we think of x as volume or mass, the measures volume and, assuming ball shaped organisms, measures surface area. So we want that substrate intake is proportional to and, ignoring metabolic costs, that g is proportional to substrate intake, hence to .

From (5.8) we deduceand next from (5.4) and (5.5)We want both g and to be nonnegative (this is an aspect that received little attention in the paper (Diekmann et al., submitted)) and so there are some restrictions. The choiceleads toand from there to the ODE systemAssuming chemostat dynamics one complements Eqs. (5.14) and (5.15) with a third equationwith the Michaelis–Menten uptake rateas a natural choice. From (5.12) can be seen to represent the influence of the resource density on the size specific cell division rate. As yet little appears to be known about this. However, given the present technology, this function can in principle be determined experimentally. Moreover, one could attempt to determine it from some model of the cell cycle in dependence on S. (Another matter is whether in reality the division rate decomposes so neatly in a product of an x-dependent and a S-dependent term.) However, our best initial guess is to take to be constant.

A catalogue of models that admit a finite dimensional state representation (a short recap of a long paper)

In Sect. 3 we already announced what is our best feat till now (Diekmann et al., submitted), a complete catalogue of all possible ODE reducible models for the case where the birth rate is treated as an output. Under this assumption we can in (3.6) drop the birth term from the backward operator as the births are already accounted for through the output map. Such models naturally allow an easy reconstruction of the full state trajectory. Additional technical assumptions were that the i-state space is an interval (possibly of infinite length) and that for all for some constant environmental condition . We then considered the problemwith w taking values in and H taking values in .

Above we have already listed the condition that (4.1) should be independent of x if so we will not return to that case (represented by model family in the paper (Diekmann et al., submitted)). For we defined two families of functions g and for which we specified w and H such that (6.1) holds. These are presented in the following catalogue.

and

Our main result is the following theorem.

Theorem 6.1

Assume that and thatThen necessarily and w can be brought into the form specified in either or by a transformation of the i-state variable and a change of basis in .

and H are such that (6.1) holds,

are linearly independent functions of .

In other words: The catalogue of solutions of (6.1) is complete.

The proof of Theorem 6.1 fills almost one half of the pages of the paper (Diekmann et al., submitted).

A few comments are in place. Firstly, it is useful to realise that a given function g may be of the form specified by even though, from the looks of it, one would not guess so. As a concrete example, considerWe challenge the reader to verify that, actually, the g given by (6.2) is of the form presented in . Readers who do not like such a challenge are invited to consult Sect. 5 of the paper (Diekmann et al., submitted).

Secondly, in the derivation of the catalogue biological constraints such as positivity of the death rate have been ignored. However, afterwards one can check whether the various families in the catalogue contain cases that satisfy positivity and other biological constraints. One of the main purposes of the catalogue is that it can be used as background information at the modelling stage.

Finally, let us recall that the Eq. (6.1) and the families are concerned with a transport-degradation model without reproduction. They extend to full population models for which (3.9) holds. In addition they are relevant for models in which all components of w represent quantities that are conserved in the reproduction process, that is, models for whichholds. To see this, simply observe that with the assumption (6.3) condition (3.7) reduces to (6.1) if one replaces by .

Physiological age and implicit scaling of time

In the developments till now we considered the requirements for a population model to be representable by an ODE whatever the initial condition. However not all applications require such a strong form of ODE representability. For example, when only some special initial conditions matter we can do with a weaker kind of ODE representability (which in analogy with similar problems in the theory of Markov chains (Kemeny and Snell 1960) we might call weak ODE representability). Note that such a special set should include also the subsequent trajectories and -limit sets. A particularly important case is when we take the latter as our special initial conditions. We shall introduce this idea by means of an example inspired by Diekmann et al. (1983) and Metz and Diekmann (1986, I.4.5, II.14). Our aim is to explain how the ideas presented there fit in with ODE reducibility as discussed above. The first step is simple. Suppose thatThen (4.1) is independent of x ifholds for some real number . We have suggestively provided with a subscript d for “dominant”. Correspondingly we assume that w is nonnegative. There exists a large literature to support this, see for instance Heijmans (1986a, b) and Bátkai et al. (2017) for general theory and the wider context. With this choice of w we find for N defined by (3.4) the ODESo far so good. The difficulty is, of course, that (7.4) provides an implicit characterisation of w and that we do not have any form of control. So, in general, this w is not suitable as a block for building a community model. Yet, as explained in detail in Diekmann (1983), Metz and Diekmann (1986, I.4.5, II.14), the variable N can be used to describe the asymptotic large time behaviour even if it fails to capture all relevant aspects of the transient phase. A crucial observation in this context is that the change of time variablefrom t to makes the problem linear, enabling us to use the theory of linear positive semigroups even though the original problem is nonlinear. The main result is that the large time behaviour of N and is described by the two-dimensional system consisting of (7.5) andThis type of asymptotic ODE reducibility was already described explicitly by Val and Metz (unpublished). They additionally considered a class of models with higher dimensional i-state space where ODE reducibility holds for an invariant subset of the population state space (while providing numerical evidence that the subset is an attractor).

Discussion

Physiologically structured population models have the advantage that one can incorporate mechanistic detail of i-level behaviour. Thus it becomes possible to investigate the relationship between mechanisms at the i-level and phenomena at the p-level (Diekmann et al. 2010; de Roos and Persson 2013). But physiologically structured population models also have disadvantages. There are only few numerical tools for their study, see Brännström et al. (2013), Carrillo et al. (2014), Breda et al. (2016), Aye and Carlsson (2017), Zhang et al. (2017), de Roos (2018) and Scarabel (2018). In this respect the situation for population models formulated in terms of ODEs is infinitely better. Moreover, there is a highly developed qualitative theory for ODEs that only partly carries over to infinite dimensional dynamical systems and when it does, answering questions concerning for instance the structure of -limit sets may be exceedingly hard. It is therefore of interest to know when a physiologically structured model can be reduced to an ODE and when it cannot.

Structured population models involve book-keeping considerations and there are multiple ways to organise these. Accordingly, there are as well multiple ways of associating a dynamical system at the population level with the model (Barril et al in preparation; Calsina et al. 2016). We favour book-keeping based on an individual’s state-at-birth and age, which leads to an abstract renewal equation for the population level birth rate (Diekmann and Gyllenberg 2007). This is particularly helpful when there is only one state-at-birth in which case the renewal equation becomes a scalar one (Diekmann et al. 2010).

In this paper we showed how the renewal equation leads, in just a few lines, to the condition (3.6) for ODE reducibility. The version (3.7) involves the per capita growth, death, and reproduction rates and as well as an valued function w of per capita outputs. After paying special attention to the cases and , in particular for fission into two equal parts, we reviewed the catalogue derived by Diekmann et al (submitted) and shown by them to be complete under certain restrictions. In addition we revealed the connection between (3.7) and the example of Diekmann (1983), Metz and Diekmann (1986, I.4.5, II.14) allowing reduction of large time behaviour to an ODE.

Originally we had some hope of being able to derive results (perhaps only partial) for the rather special, but also important, caseof only one possible state-at-birth . In fact we got nowhere. So we like to present this as an open problem.

F_2: Physiological age
k	k≥2
Parameters	γ_0:Ω→R,μ_0:E→R,v_0:Ω→R,v_1:E→R
	Λ∈R^k×k and w(x_b)∈R^k
	ζ(x):=∫x_bxdyv_0(y)
g	g(x,E)=v_0(x)v_1(E)
μ	μ(x,E)=γ_0(x)g(x,E)+μ_0(E)
w	w(x)=exp∫x_bxγ_0(y)dyexpζ(x)Λw(x_b)
H	H(E)=v_1(E)Λ-μ_0(E)I

F_3: Generalised von Bertalanffy growth
k	k≥2
Parameters	γ_0:Ω→R,μ_0:E→R,v_0:Ω→R,v_j:E→R,j=1,2,3.
	ζ(x):=∫x_bxdyv_0(y)
g	g(x,E)=v_0(x)v_1(E)+v_2(E)ζ(x)+v_3(E)ζ(x)^2
μ	μ(x,E)=γ_0(x)g(x,E)+μ_0(E)+(k-1)v_3(E)ζ(x)
w	w_j(x)=exp∫x_bxγ_0(y)dyζ(x)^j-1,j=1,2,…,k.
H	H(E)=H_0(E)-μ_0(E)I
	H_0=0-(k-1)v_30⋯0v_1v_2-(k-2)v_3002v_12v_20⋮03v_10⋱0⋯⋯(k-1)v_1(k-1)v_2