## Abstract

The quick answer to the title question is: by bookkeeping; introduce as p(opulation)-state a measure telling how the individuals are distributed over their common i(ndividual)-state space, and track how the various i-processes change this measure. Unfortunately, this answer leads to a mathematical theory that is technically complicated as well as immature. Alternatively, one may describe a population in terms of the history of the population birth rate together with the history of any environmental variables affecting i-state changes, reproduction and survival. Thus, a population model leads to delay equations. This delay formulation corresponds to a restriction of the p-dynamics to a forward invariant attracting set, so that no information is lost that is relevant for long-term dynamics. For such equations there exists a well-developed theory. In particular, numerical bifurcation tools work essentially the same as for ordinary differential equations. However, the available tools still need considerable adaptation before they can be practically applied to the dynamic energy budget (DEB) model. For the time being we recommend simplifying the i-dynamics before embarking on a systematic mathematical exploration of the associated p-behaviour. The long-term aim is to extend the tools, with the DEB model as a relevant goal post.

## 1. Introduction

Within the framework of physiologically structured population models (PSPM) one can, in principle, incorporate a lot of mechanistic details about physiological processes at the i-level (i for individual), such as that found in the dynamic energy budget (DEB; Kooijman 2009*b*) models, which form the subject of this theme issue, and yet arrive at a consistent and complete deterministic bookkeeping scheme for a sufficiently large population (Metz & Diekmann 1986). The aim of such exercises is to deduce how population phenomena relate to the mechanisms of feeding, metabolism, maintenance, growth, reproduction, starvation-induced death, and so on. But in order to carry out this deduction, one needs a mathematical framework that provides the tools.

The conceptual part of the mathematical framework is easy: once the notion of ‘state’ has been introduced at the i-level (in the DEB model structural mass, non-allocated reserve mass, and so on), it can be lifted to the p-level (p for population) by declaring the p-state to be the measure *m* on the i-state-space *Ω*, such that *m*(*Ω*) is the total population size and for every measurable subset *ω* of *Ω*, the fraction of the population with i-state in *ω* is *m*(*ω*)/*m*(*Ω*)—in other words, *m* describes the population size and distribution—and cataloguing how the various i-processes contribute to the change of *m*.

The next task is to check that the resulting bookkeeping scheme constructively defines a dynamical system, that is, unambiguously defines the measure *m _{t}*

_{+s}, describing the population at time

*t*+

*s*, when the measure

*m*, describing the population at time

_{s}*s*, is specified. While carrying out this task, one is confronted with a surprising amount of technical difficulties (Thieme 1988; Diekmann

*et al*. 2001; Diekmann & Getto 2005), many of which reflect subtle modelling issues (for example, if the rate of channelling energy to reproduction changes abruptly when an individual passes a certain size, as it seems to do for

*Daphnia*, what does an individual do when it happens to stop growing exactly when reaching this size (Thieme 1988)?). Addressing these issues has some spin-off in terms of increased biological/modelling insight. Yet, given the fact that the construction of a dynamical system is a very preliminary contribution to unravelling the p-behaviour, it is somewhat depressing that so much hard work is needed for what ought, one feels, to be a simple task.

Things turn really bad if one attempts to develop stability and bifurcation theorems and tools. The chief difficulty is a severe lack of smoothness if the i-growth rate is affected by, for instance, competition for food. (Smoothness is here meant in a very abstract sense: the map from food availability as a function of time to future population states is in general not even differentiable thanks to the untoward geometry of the space of measures.) Even if one can derive a characteristic equation by formally linearizing around a steady state and then looking for exponential solutions (e.g. de Roos *et al*. 1990), a rigorous proof that the positions of the roots of this characteristic equation do indeed govern the (in)stability of the steady state could not be found.

Sometimes one can work around a problem that one cannot solve (perhaps because it *is* unsolvable). The aim of this paper is to demonstrate that, fortunately, such is the case for PSPM. The key point is to reflect on the notion of ‘state’ and to see whether it can, possibly after a suitable and justified restriction, be represented in a different way.

The ‘state’ of a system is, by definition, the information about the past that is relevant for predicting the future. A population ‘experiment’ is started by specifying the states of the individuals in the ‘inoculum’. If food supply and predation pressure (and possibly other characteristics of the world in which these individuals live) are given functions of time, maturation, survival and reproduction can be ‘computed’. Once they are born, the same can be done for the individuals that arise from reproduction. In fact, however, quantities like food supply and predation pressure are not given functions of time, as they are influenced by the focal population. This we call feedback by way of environmental interaction variables. The prediction of the future hence involves the solution of a coupled system of equations for the p-birth rate (taking values in the space of conceivable states-at-birth, which usually is much smaller than the full i-state space), and the environmental interaction variables, all as functions of time. Once these are known, one can compute how the measure describing the population size and composition depends on time (Thieme 1988; Diekmann *et al*. 2001).

It follows that, at the p-level, the information about the past that is relevant for predicting the future, is the *history* of the p-birth rate and the environmental interaction variables. In this formulation, we deliberately ignore the individuals from the inoculum that are still alive. ‘Deliberately’, since their contribution fades with time. So we make the restriction that we only consider those measures that can be constructed from the history of the p-birth rate and the history of the environmental interaction variables. This restriction is justified, in the sense that these measures constitute a forward invariant subset of the p-state space that attracts all orbits (if individuals have a uniform upper bound to their age *a*_{max}, in state-at-birth and feasible environmental conditions, then the dynamical system actually maps any initial condition into the invariant set for time bigger than *a*_{max}). This restriction is suitable, in the sense that the lack of smoothness is eliminated and therefore a rich stability and bifurcation theory can be developed with some, but not too much, effort (cf. Diekmann *et al*. 1995, 2007; Diekmann & Gyllenberg submitted). (There is some wishful thinking in this formulation: if the behaviour of individuals changes abruptly upon passing a certain size, one needs the more involved theory concerning state-dependent-delay equations; cf. Hartung *et al*. 2006.)

Once we adopt the aforementioned restriction, we may actually shift the measures to the background and put, at the p-level, the spotlight on the p-birth rate and the interaction variables. This we call the delay equation formulation of PSPM.

The plan of this paper is as follows. In §2, we illustrate the delay equation formulation by presenting the example of a population structured according to a one-dimensional i-state, for definiteness to be called size, feeding on an unstructured resource, while referring to Diekmann *et al*. (2010) for details and results. In §3, we briefly sketch how to build numerical bifurcation tools for the analysis of p-models formulated as delay equations, this time referring to de Roos *et al*. (2010) for a full exposition. In §4, we discuss the potential for applying the ideas sketched in the previous section to the DEB model, and in the final section we discuss how the new framework relates to earlier ones, as well as our expectations for the future.

## 2. The delay equation formulation

In this section, we formulate a model of the interaction between an unstructured resource and a consumer population structured according to a one-dimensional physiological variable that we shall refer to as ‘size’. We assume that all consumers are born with the same size *ξ*_{b} and that their growth is deterministic according to the differential equation
2.1
where *ξ* denotes size, *a* age, *X* resource concentration and *t*(*a*) the time at which the focal individual has age *a*. We assume that the survival probability of an individual decreases according to
2.2
and that newborns are produced clonally at a rate
2.3

The energy and materials needed for maintenance, growth and reproduction are derived from the ingestion of the resource, which proceeds at a rate
2.4
but for the time being we leave the relationship between, on the one hand *g*, respectively, *β* and, on the other hand *γ*, unspecified. Finally, we assume that in the absence of consumers the resource evolves according to the differential equation
2.5

Essentially, the model is now specified, the rest is bookkeeping. To do the bookkeeping in an efficient way, we need to introduce some notation. As in the theory of delay equations, we use the symbol *X _{t}* to denote the history of the food concentration relative to

*t*, i.e. the function 2.6 (In case of a maximal age

*a*

_{max}, we can restrict to

*σ*≥ −

*a*

_{max}, but otherwise we need to allow

*σ*to take any value less than or equal to zero.)

Now consider an individual that has age *a* at the current time *t*. Denote the size of this individual at age *τ*, with 0 ≤ *τ* ≤ *a* by *ξ*(*τ*) = *ξ*(*τ*;*a*,*X*_{t}) (the last two variables are suppressed in the notation whenever that helps to keep formulae readable). It can be computed from
2.7
(Age is counted here from the moment of birth, most often in the form of an egg being laid; what matters mathematically is that the link between mother and young is severed.) The size at the current time is then given by
2.8

The fraction of individuals born at time *t* − *a* that are still alive at time *t*, is given by the survival probability . Let be the survival probability up till time *t* − *a* + *τ*. Then, by definition,
2.9
with computed from
2.10

Let *b*(*t*) denote the consumer population birth rate at time *t*. Then the assumptions above imply that
2.11
The upper equation is the time-dependent analogue of Lotka's equation, with the *pro capite* birth rate and survival continuously updated in dependence on the experienced environmental history. The integral in the lower equation adds up the rates of resource consumption by individuals of different ages.

As written, equation (2.11) refers to consumers and a resource that have been interacting since the dawn of time. When looking for steady states, periodic solutions and so on, this is indeed the right perspective. But alternatively, one can require the system of equations (2.11) to hold only for positive time and provide initial data in the form of the history of both *b* and *X* at time zero, the first as a non-negative locally integrable function and the second as a non-negative continuous function. In the mathematical theory concerning equation (2.11), for which we refer to Diekmann *et al*. 2007 (in particular, §3), both points of view play a role.

Note that the right-hand side of equation (2.11) is linear in *b*, reflecting that the environmental interaction variables are chosen such that all dependence is by way of feedback via these variables (i.e. if one considers *X* as prescribed, then the consumers are independent of one another).

We now consider constant solutions of equation (2.11), i.e. steady population states. If food is kept at the constant concentration , then the basic reproduction number *R*_{0} of the consumers is well defined and given by
2.12
(where we use the same symbol to denote the value and the constant function taking that value). The equation
2.13
determines the constant concentrations that lead to a steady consumer population. As a rule, any solution of equation (2.13) is unique, simply since *R*_{0}(0) < 1, *R*_{0}(*∞*) > 1 and *R*_{0} is a monotone function of . (Of course, other scenarios are possible, for example food that is toxic at high concentrations.)

Once is determined, we have explicitly
2.14
for the consumer population birth rate that keeps the resource at the steady level . (The steady states of more general PSPMs can be determined by following essentially the same steps; see Diekmann *et al*. (2003)).

To linearize equation (2.11) around the steady state is, in principle, straightforward, but very laborious, so we chose not to give the details; for these see Diekmann *et al*. (2010). One obtains a system of the form
2.15
The corresponding characteristic equation is
2.16
or
2.17
where the hat denotes Laplace transform, that is,
2.18
Equation (2.18) is obtained by looking for exponential solutions of equation (2.15).

The hard work consists of

— computing the ingredients

*c*and*k*(*a*) of equation (2.15) from*g*,*μ*,*β*,*γ*and*h*,— analysing equation (2.17) in order to determine conditions on

*g*,*μ*,*β*,*γ*and*h*, which guarantee that all roots are in the left half of the complex plane and complementary conditions that guarantee that at least one pair of roots lies in the right half of the complex plane.

The results derived in Diekmann *et al*. (2007) imply that one can draw conclusions about the dynamical behaviour of solutions of equation (2.11) from this kind of information about the position in the complex plane of the roots of the characteristic equation. In Diekmann *et al*. (2010), the first point is elaborated in quite some detail and analytical as well as numerical results concerning the second point are presented. In general, however, one needs computer assistance when doing the hard work. So, in the next section we focus on some computational tools that have been developed in order to carry out an analysis of equation (2.11) and its linearization in equation (2.15).

## 3. The existence and stability boundaries

In order for equation (2.14) to be biologically meaningful, we need that , where is such that equation (2.13) holds. The borderline is determined by the system of two equations 3.1 in the single unknown . Of course, as already expressed by the word ‘borderline’, there is generically no solution. But when we ‘free’ one parameter (i.e. consider this parameter as unknown; its choice to be determined by a suitable combination of convenience and biological relevance) we might be able to solve equation (3.1). Yet, we recommend to free two parameters and to consider equation (3.1) as two equations in three unknowns such that, generically, there is a curve of solutions. (We have two arguments for this recommendation: (i) curves can numerically be efficiently computed by means of continuation methods; and (ii) humans are very well equipped to absorb information that comes in the form of a two-dimensional picture.) The projection of this curve on the two-dimensional parameter space is called the existence boundary, since it separates the region for which equation (2.13) has a biologically meaningful solution from the region where it has not (and where the consumer is doomed to go extinct).

Typical examples of functions *h* are

— chemostat resource dynamics:

*h*(*X*) =*D*(*X*_{0}−*X*),— logistic resource dynamics:

*h*(*X*) =*rX*(1 −*X*/*X*_{0}).

A typical choice of the two parameters is *X*_{0} and a uniform consumer death rate *μ*_{0}. The equation then leads to , and the existence boundary is determined by the single equation
3.2
To trace numerically the curve defined by equation (3.2) (for instance, by a predictor–corrector method), we need, first of all, to be able to evaluate *R*_{0} for given *X*_{0} and *μ*_{0}. This can be done by solving the system of ordinary differential equations (ODEs; e.g. de Roos 2008)
3.3
with *R*_{0} = *B*(∞). Of course, we do not want to integrate equation (3.3) forever. A stopping criterion can often be based on some variant of the following result: if we integrate only to some *a*_{max}, then, if *g*(*ξ*; *X*_{0}), *μ*_{1}(*ξ*; *X*_{0}) and *β*(*ξ*; *X*_{0}) are monotone beyond *ξ*(*a*_{max}), *g*(*ξ*; *X*_{0}) decreasing, *μ*_{1}(*ξ*; *X*_{0}) and *β*(*ξ*; *X*_{0}) non-decreasing and there exists a , such that , then
3.4
Both equation (3.4) and the fact that *g*, *μ*, *β* and *γ* may have jump discontinuities at *ξ* = *ξ*_{A}, the size at which juveniles turn adult (meaning that they start to reproduce), imply that we need to incorporate in the ODE integration routine criteria that test whether a certain variable is still below a threshold.

Referring to de Roos *et al*. (2010) for more details, we conclude that one can define the left-hand side of equation (3.2) in terms of a procedure that solves the ODE system (3.3) and next compute the existence boundary in the (*X*_{0}, *μ*_{0}) plane by solving equation (3.2) numerically with a standard continuation method.

At the existence boundary, a transcritical bifurcation takes place. In the absence of Allee effects, the bifurcation is supercritical, which means that within the existence area close to the boundary, a steady state exists with *b* as defined by equation (2.14) being small. According to the Principle of the Exchange of Stability (see Boldin 2006, and the references therein) this steady state is stable. The steady state may lose stability further on in the existence area. The curve that separates the stability area from the instability area is called the stability boundary. It can be computed in much the same way as the existence boundary, but this computation is a lot more complicated, as now the linearization (3.1) is involved as well.

In the absence of fold bifurcations (as is the case when the solution of equation (2.13) is a single-valued function of the parameters), the stability boundary is characterized by the characteristic equation (3.3) having a purely imaginary root *λ* = *i**ω*, *ω* ≠ 0. The characteristic equation is a single complex equation, hence counts for two real equations. So, together with equation (2.13), this amounts to three equations in the four unknowns , *ω*, *X*_{0} and *μ*_{0}. The strategy to find a solution curve is as before first to find one point on the curve and next use a continuation algorithm to compute the entire curve. But now we need to evaluate, for given , *ω*, *X*_{0} and *μ*_{0}, not only *R*_{0} but also the ingredients of the characteristic equation. This means that we have to integrate equation (3.3) but also
3.5
to find the cumulative food consumption, which is needed to evaluate the right-hand side of equation (2.14), and the linearized versions of these ODE with forcing terms that involve cos(*ω**a*) and sin(*ω**a*). Details and pseudocode are presented in de Roos *et al*. (2010).

## 4. Extensions and applications to the deb model

The mathematical theory referred to in §1 is in principle very general. However, that generality only holds good at an abstract level. The more practical the results one aims at, the more restrictions have to be imposed in order to get results. At the more down to earth end, we have so far worked out the concrete details only for the small subset of simple cases discussed in §§2 and 3. So we can presently do little more than indicate the difficulties and their potential solutions that we envision for handling DEB-style models in the rigorous mathematical fashion that we are aiming at. In doing so, we shall use a mix of our earlier notation and that of Sousa *et al*. (2010); we only explain conventions when there is a clash or when we have to introduce additional symbols.

Difficulties in applying a mathematical framework can be of very different kinds. First there is the problem of whether a concrete description of individual behaviour can be represented at all within an envisaged framework for population modelling. The DEB model no doubt fits in the very general physiologically structured population framework. Then there is the problem of whether a model falls within a class for which at least the mathematical existence and uniqueness problem for the heuristically derived population equations can be solved. The description of DEB models certainly appears compatible with the existence and uniqueness results in Diekmann *et al*. (2001), and in principle also with the framework described in §§2 and 3. However, they clearly do not fit the example that formed the centrepiece of those sections. As the results for that example only exemplify the more general theorems in Diekmann *et al*. (2007) and Diekmann & Gyllenberg (2008, submitted), the first question then is to what extent DEB models fit those abstract theorems, and if not, whether they, or the theorems, can be tweaked to make them fit. A second question is then how far the abstract results can in principle be implemented in the form of manageable calculations. And a final question is what are the best strategies to implement those calculations. We will not consider the questions in this order but rather use various technical specifics of DEB models as our guideline.

The natural strategy for answering the sort of questions indicated above starts with writing the DEB model presented in Sousa *et al*. (2010) in a form similar to the one used in §2, in order to see if it indeed satisfies the requirements of the theorems that have been proven this far. Below we shall indicate how this should be done. It is more practical, however, first to introduce a simplification meant to remove an aspect of the DEB models that from a mathematician's point of view is technically somewhat worrisome. The full DEB model is rigidly deterministic at the level of the individual. In the DEB-model proper, young are produced from the reserves set aside for reproduction *M*_{ER} by, when *M*_{ER} reaches a threshold, converting those reserves into young while resetting *M*_{ER} to zero. Whenever the i-level model is consistently specified, it is, of course, always possible to do individual-based simulations. However, our goal is to analyse the deterministic models that result as large number limits of such stochastic descriptions. In general, the mathematical tractability of these deterministic models to a large extent depends on the fact that over time too sharp bumps in the population state get smoothed out. Not only that, for a given course of the environment, *X*(*t*), *t* > 0, the population states over time will look more and more like each other, but for a difference in the total population sizes. (There are exceptions, which by this very fact are of considerable mathematical but little biological interest.) Having the young produced deterministically tends to thwart this smoothing property to an extent that the main present proof techniques fail. This does not mean that there is something mathematically wrong with the full DEB model. Only that the techniques are not yet up to this complexity. Experience teaches that often such problems have to do with the existence of badly behaving exceptional cases that for all practical purposes can be neglected. ‘In reality’ all model ingredients tend to be subject to minor chance fluctuations that may be expected to remove the mathematical anomalies if incorporated in the model structure. However, doing so leads to models of severely increased complexity. In order to keep a somewhat tractable framework, we in the past fudged the effect of stochasticity in the individual state transitions by having the process of giving birth represented as a continuous rate *β* (cf. Metz & Diekmann 1986, Remark I.3.2.1; see their Section III.6.3 and Heijmans & Metz 1989 for one possible justification, in the form of a limit argument starting from a more realistic specification). This is also what we did in §§2 and 3 and what we will do below: we shall assume that births are produced at a rate that is proportional to the rate at which energy is channelled towards reproductive activities.

An additional advantage of assuming that the birth process is continuous, is that this way one state variable, *M*_{ER}, is removed, which also is the one that does not fit too happily in the procedure for calculating a characteristic equation outlined in §2 owing to its jumping behaviour. In principle, this difficulty can be overcome through calculations like those in Metz & van Batenburg (1985, Section 4.2) and Kooi & Kooijman (1999), but the characteristic equation will become rather complicated. This practical aspect comes on top of the mathematical problem that, owing to the difficulties mentioned earlier, one cannot be fully sure that information from the characteristic equation indeed has the implications that we normally attach to it (although at least one of us dares to bet on it).

With the assumption about the birth process that we just made, the DEB model as described in Sousa *et al*. (2010) has a seven-dimensional i-state, which we choose to represent by , where the state variables and represent the maxima of *M*_{V} and *M*_{H} along their trajectories till the present age.

The rates of change of the first five state variables are given by equations (1, 8, 9, 12, 13, 14) in Sousa *et al*. (2010). The latter two variables satisfy
4.1

The death rate equals
4.2
*θ*_{H} and *θ*_{V} parameters, where the inequality conditions exclude rejuvenation and shrinking due to starvation (van Leeuwen *et al.* 2002; Kooijman 2009*b*).

The *pro capite* birth rate *β* can be found in Sousa *et al*. (2010) in the section on maturation and reproduction under the name . At the birth of an egg, , where is such that *M*_{E}/*M*_{V} at the start of the juvenile stage () equals that of the mother at egg formation (details in Kooijman 2009*a*).

The numerical integration of the corresponding population equations, e.g. with the Escalator Boxcar Train method of de Roos (1988, 1997), de Roos & Metz (1991) and de Roos *et al*. (1992) (see also http://staff.science.uva.nl/~aroos/), runs into a number of practical problems that all have to do with the determination of the birth state. The maternal effect from the standard DEB model can cause young born at the same time to be distributed over a continuum of birth states. We expect, however, that in most applications the birth states of those young will be sufficiently close to each other that lumping them in but a few classes, to be represented by their class mean, will not cause much of a discretization error. An additional complication is caused by the lack of an explicit expression for , which hence has to be determined by numerically solving an equation (an efficient numerical scheme may be found in Kooijman 2009*a*). However, none of these problems is fundamental.

The previous developments show that it is possible to express the DEB model in a form compatible with the ideas from §§2 and 3. Another matter is that the DEB equations are far more complicated than the simple equations to which we restricted ourselves in §2, so that one cannot immediately apply the concrete recipe given there. Below we shall systematically discuss any technical problems that this complication may give.

The fact that in the DEB model the i-state is higher dimensional in principle does not greatly complicate the calculation of the kernels *k _{ij}* from §2. One only has to interpret the various scalar expressions from which the

*k*are calculated as vectors and matrices. (The last two state variables can be dropped from consideration since close to equilibrium the inequality conditions in equation (4.2) are never violated.)

_{ij}A greater problem is the multiplicity of birth states. These even do not come in a finite number but in a continuum. Rigorously extending the proofs in Diekmann *et al*. (2007) to this case brings considerable technical difficulties with it. At a more practical level, the local stability of systems with infinitely many birth states can in general no longer be determined from the analysis of a characteristic equation. For the case of finitely many birth states, it is possible in principle to extend the calculations of Diekmann and co-workers by further invoking the vector–matrix formalism. All that results is that the matrix in equation (2.16) becomes correspondingly larger. Where equation (2.16) contains a 2 × 2 matrix, with *n* birth states, that matrix will become an (*n* + 1) × (*n* + 1) one. (When, as in the DEB model, the attribution of birth states has to be solved from an equation, one has to differentiate through that equation and solve the resulting linear equations for the required derivatives.)

A similar extension holds if the number of environmental variables, be they resources, toxicants or predation pressures, is larger than one. Assuming for the time being that all these variables satisfy simple differential equations, with *m* environmental variables the matrix will become of size (*n* + *m*) × (*n* + *m*). Of course, with any increase of *n* or *m*, the calculations will become more forbidding and will ultimately fail when either number becomes infinite. (An example where the number of environmental variables is infinite are models where a resource has a continuous size distribution, with the speed of ingestion depending on the form of this distribution in a non-trivial manner, as is the case in various published size structured predator–prey models.)

It may be expected that even the cases with a continuum of birth states and/or environmental variables may ultimately be handled by some discretization, such as is anyway done in numerical calculations. However, theorems that the stability bounds found by such approximate calculations approximate the ones of the original model are still pending.

In short, we are still a long way from tackling the full complexity of the DEB model in a rigorous fashion outlined in §§2 and 3. Given their well-specified mathematical structure and their great and still manageable (e.g. through individual-based simulation) biological realism, DEB-style models provide an interesting pointer, as well as set of test cases, for how to proceed.

## 5. Discussion

We have in the past (Metz & Diekmann 1986) promulgated representing structured population models by partial (integro-)differential equations (PDEs) for the density of the population over its i-state space, like in Kooi & van der Meer (2010) for example. In hindsight, in a strict mathematical sense a lot turned out to be wrong with this formulation, as it can be interpreted rigorously only for a very special subset of the biological systems for which it was envisioned as representational. However, in the cases we have scrutinized so far, the validity of the biological conclusions arrived at this way is not affected by these mathematical difficulties (see below).

A first problem with the PDE framework is that in theory it is possible to start with a cohort of equal individuals so that there is no population density to begin with. The same problem may occur in the probably more familiar problem of diffusive movement in space. However, there this initial singularity disappears immediately after *t* = 0. In the case of structured populations with deterministic i-state movement, these point masses stay intact till the ancestral population has died out. So, if we are interested in asymptotic calculations only, there is not much to worry about. A bigger problem is that for many models the mechanism itself creates distributions over the i-state space that do not admit a density, like when all individuals are born equal and move in exactly the same manner through a *k*-dimensional i-state space, *k* > 1. In that case all individuals born after *t* = 0 are concentrated on a one-dimensional curve in a *k*-dimensional space. Perhaps, these technical anomalies could still be handled by using some sort of ‘weak solution’ concept. However, we found it easier to start from a representation of the population state as a measure, which is what it ultimately is anyway!

On top of this then comes the problem of the dependence of the differential equation describing the i-state movement on the current state of the system, consisting of the population and, say, its food. We already discussed in the introduction how this biologically reasonable assumption gets us in deep waters mathematically. In Diekmann *et al*. (2001), we therefore proved the existence and uniqueness by interpreting the population equations as an input–output operator from environment to environment and solving the fixed-point problem that appears when we connect input and output. (Solving here means constructing the solution by means of abstract mathematical operations, like integrations and taking limits, which is different from arriving at an explicit expression or constructing an efficient numerical algorithm.) However, that is as far as we got.

We believe that we have now found a way to the shallows again by moving from a dynamics on a space of measures to one on a space of histories of the birth rate and environment. The latter spaces have a more amenable structure, both since they represent some smoothed representation of the population state and since they look only at the subset of the population states that remain after the effects of all too bumpy initial conditions have worn off. Yet, this space is still sufficiently large that a meaningful stability theory can be built on it.

The fact that the PDE formalism cannot be considered mathematically fully sound in practice detracts but little from the papers seemingly based on that formalism, for these actually never use the formalism in the form in which they present it. For example, they use numerical techniques like the Escalator Boxcar Train for the following population trajectories over time (de Roos 1988, 1997; de Roos & Metz 1991; de Roos *et al*. 1992), which is eminently compatible with the delay equation formalism (but for the representation of the inoculum, which however enters in the main scheme in a compatible way, as an influence on the resource dynamics and a birth stream), or they formally calculate a characteristic equation, which can be seen to be the same as they would have got from the delay equation formalism, except that they move in a different order through the various calculation steps. (More in particular, a formal integration along characteristics, which one can do either in the PDE or after a few of the standard steps for arriving at a characteristic equation, to wit inserting exponential test solutions into the equation and separating terms, leads precisely to the ingredients of the delay equation formalism, as this integration is no more than the following of individuals over their lives from which the latter formalism is derived in a direct manner. However, beware: handling discontinuities in the coefficient functions can be rather tricky in the PDE formalism. Here, the delay equation formalism proves its mettle also for practical calculations.)

This ends the exposition of the technical mathematical problems and ways to overcome them. For a finale, let us take one last look at how the various model structures connect to real ecology. We seemingly failed to take account of the spatial extent of populations. In principle, this can be remedied by taking space as an extra i-state variable. In the delay equation formalism, space would come in as an additional component of the birth state. So in principle not much changes. However, practically there are two problems. One is that in the general case it is no longer possible to derive a characteristic equation. Only for particular simple shapes of the spatial domain (amenable to a description with a separable coordinate system), it may be possible to calculate a so-called dispersion relation, which in essence fulfils the same role. The second difficulty is that movement in space is generally modelled as random. In the simplest case, where we have only a finite number of locations each supporting a well-mixed population, calculations for random movements, similar to those for the deterministic movements of the i-state that we considered so far, probably are rather painless, but at the time of writing even that work still has to be done.

A further lack of realism is that we only considered a single species in isolation, but for its dependency on a dynamical food source. The general existence and uniqueness results of Diekmann *et al*. (2001) immediately apply to the multi-species case. We already indicated how the calculation procedures extend in the case where the other species can be represented as scalars, similar to the resource in our ‘single’ species model. For more than one structured species, we may expect the calculations to go along similar lines, provided the number of interaction variables can be kept finite. However, the work needed to get anywhere useful increases very fast with every increase in model complexity.

In summary, in principle, the ideas proposed in §§2 and 3 allow a fair amount of extension. In practice, it is for the time being probably better first to simplify the models considerably before even asking your local mathematician to work on them, as mathematicians also have to learn their trade by increasing the complexity of the problems they work on in small steps.

## Footnotes

One contribution of 14 to a Theme Issue ‘Developments in dynamic energy budget theory and its applications’.

- © 2010 The Royal Society