Atto-Foxes and Other Minutiae

This paper addresses the problem of extinction in continuous models of population dynamics associated with small numbers of individuals. We begin with an extended discussion of extinction in the particular case of a stochastic logistic model, and how it relates to the corresponding continuous model. Two examples of ‘small number dynamics’ are then considered. The first is what Mollison calls the ‘atto-fox’ problem (in a model of fox rabies), referring to the problematic theoretical occurrence of a predicted rabid fox density of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-18}$$\end{document}10-18 (atto-) per square kilometre. The second is how the production of large numbers of eggs by an individual can reliably lead to the eventual survival of a handful of adults, as it would seem that extinction then becomes a likely possibility. We describe the occurrence of the atto-fox problem in other contexts, such as the microbial ‘yocto-cell’ problem, and we suggest that the modelling resolution is to allow for the existence of a reservoir for the extinctively challenged individuals. This is functionally similar to the concept of a ‘refuge’ in predator–prey systems and represents a state for the individuals in which they are immune from destruction. For what I call the ‘frogspawn’ problem, where only a few individuals survive to adulthood from a large number of eggs, we provide a simple explanation based on a Holling type 3 response and elaborate it by means of a suitable nonlinear age-structured model.


Introduction
In a recent fascinating conversation with the ecologist Yvonne Buckley, of Trinity College, Dublin, we touched on a number of issues concerning population dynamics which have puzzled me for some time. In this paper I want to draw these conundrums together and offer some palliative solutions.
The central theme is the use of continuous population models in circumstances where the population levels become extremely low. A continuous population model relies on the assumption that the population size varies continuously in time. This requires, for two reasons, that the population size be large. Firstly, because the actual discrete changes in integer numbers need to be viewed as infinitesimal changes, and secondly, because actual finite time gestation periods can only be interpreted as continuous in time if the large population can be taken as a realisation of an evolving probability density. Continuous population models are in particular subject to criticism when they indicate very low population levels, and in this circumstance discrete and/or stochastic models may be preferable (Durrett and Levin 1994).
The issue is simply illustrated by a linear birth-death process (Bartlett 1960) in which, if the population is of size n, individuals have probability λ n δt of giving birth in a time interval δt and equivalent probability of dying of μ n δt. If p n (t) is the probability that the population is of size n at time t, then it is straightforward to show thatṗ n = (n − 1)λ n−1 p n−1 − n(λ n + μ n ) p n + (n + 1)μ n+1 p n+1 ; (1.1) this applies for n ≥ 0, providing we define p −1 = 0. Note that p 0 is the probability of extinction at time t. If λ n = λ and μ n = μ are constant, the equation is easily solved with a generating function (1.4) The mean of the population is and it is simple to show directly from (1.3) that N satisfies the mean field equatioṅ (1.6) The probability of extinction, p 0 = G(0, t), is given by ( 1.7) and we see that If the population is growing (λ > μ) and the initial population size is of reasonable size (m 1), then the likelihood of extinction is negligible; on the other hand a decaying population (λ < μ) will eventually become extinct. Essentially the same conclusion is true if λ and μ are functions of t.

Nonlinear Stochastic Models
The problem with such discrete models is that their extension to nonlinear processes becomes less tractable. As the simplest example, suppose that the birth rate λ is constant, but that the specific death rate is μ n = λn K . In this case we would expect that the mean of the population would satisfy the Verhulst (1845) logistic equation, with carrying capacity K . There is a large literature dealing with such problems; see for example Bartlett et al. (1960), Nåsell (2001), Ovaskainen and Meerson (2010) and Doering et al. (2005). The paper by Ovaskainen and Meerson, in particular, contains many further references. We summarise some of the results here, though perhaps in a slightly different guise.
The equation (1.1) takes the forṁ 9) and the generating function defined in (1.2) now satisfies where we have scaled time by writing τ = λt. As is commonly done, the subscripts τ and s in (1.10) denote partial derivatives (but the subscripts n, etc. in (1.9) are indices).

Generating Function
The carrying capacity K is an integer, and if it is not large, we would expect extinction to occur fairly rapidly. In fact, eventual extinction is certain for any finite K , in the sense that G → 1 as τ → ∞ (thus p 0 (∞) = 1). On the face of it, this seems to indicate that a continuum model is doomed. If we consider (1.9) for n = 0, 1, . . ., it is not difficult to show, since p 0 is bounded above and thusṗ 0 → 0 as τ → ∞, that p n → 0 for all n ≥ 1. The two possibilities are that the mean of the population grows unboundedly, or that extinction occurs. The effect of the diffusion term in (1.10) appears to imply the second of these conclusions. If we write G = 1 + g (so g also satisfies (1.10), but its steady state is g = 0), then we find (1.11) and thus indeed g → 0 (since g = 0 at s = 1). The rate of approach can be estimated by writing g = φ(s)e σ τ , and then conversion of (1.10) to Sturm-Liouville form where the admissible functions φ are piecewise smooth functions with φ(1) = 0. A simple estimate (which also betrays the boundary layer structure of the eigenfunctions) is to take (1.14) whence we find for large K that −σ < ∼ 1.54K e −K (1.54 is the minimum of e λ − 1 λ 2 ). From this we see that for large K , extinction takes an exponentially large time to occur. This is similar to the Ehrenfest urn problem and is not a matter for concern when K is large, but for O(1) values of K , extinction will occur on times τ ∼ O(1). In practice, the distribution approaches a quasi-steady state (Bartlett et al. 1960), which may be determined as follows. In a steady state, (1.9) implies p n = (n − 1)K n 2 p n−1 , n ≥ 1, (1.15) and obviously p 0 = 1, p n = 0 for n > 1. The quasi-steady-state assumption is that p n for n ≥ 1 approaches a (quasi-)steady state while p 0 < 1; the idea is that p 0 varies on a slow time scale. If this is the case then we can solve (1.15) to obtain where A is slowly varying. Applying G = 1 at s = 1, this gives ( 1.17) we can then use this to calculatė using Laplace's method, which confirms the quasi-steady-state hypothesis and is also consistent with the earlier estimate of decay rate. The next question is to provide this quasi-stationary profile. This can be obtained from (1.16) using Stirling's approximation, but it is more illuminating to use a continuum approximation for p n in order to show that it is obtained on a timescale of τ ∼ O(1).

Continuum Approximation
Keeping K large (when we therefore might expect the continuous model to apply), we revert to (1.1), and then writing (1.1) is a discrete approximation to which can be solved using the method of characteristics, for example with an initial condition p = δ(ξ − ξ 0 ), where K ξ 0 is the initial population size. It can then be shown that the continuous Verhulst model for the mean population is regained. A nice way to demonstrate this uses the fact that with a delta function as initial condition, the solution must in fact be (1.21) Substituting this in to (1.20), we find, using the language of generalised functions, now we multiply by (ξ − N ) and use the fact that xδ(x) = 0 and thus xδ ( on integrating over 0 < ξ < ∞. An extension to this for large K is to use the next term in the expansion of (1.1), which leads to the Fokker-Planck equation 24) and to solve this, we write and if we choose N to satisfy (1.23), then (1.26) As τ → ∞, N → 1, and (1.26) has a quasi-steady solution There is a caveat to this result. This is because the approximation in (1.27) is invalid for η ∼ √ K . To deal with this, we revert to p and ξ , and since the far-field expression in (1.27) is (1.30) Using the fact that (1.31) The initial condition we choose for (1.31) should correspond approximately to an initial delta function, which we take to be centred at the steady value ξ = 1. Note that an arbitrary additive constant for ψ can be absorbed into A. To represent the initial condition, we will consider the family of functions where the limit of a → ∞ corresponds to a delta function. which can be solved by writing it in characteristic form using Charpit's equations. The initial condition (1.32) can be written in parametric form as and Charpit's equations reduce to where the overdots refer to differentiation with respect to τ . The expression for P comes from solving the quadratic equation for e P given by F = 0. The upper or lower signs in (1.35) are chosen so thaṫ This function is plotted in Fig. 1, along with Q 0 (σ ). Evidently the upper sign is chosen for 0 < ∼ ξ < 1 and the lower one for ξ > 1. Thus for small τ , the characteristics move to the left for ξ < 1 and to the right for ξ > 1. This will remain true unlessξ = 0, which occurs when Q = Q ± (ξ ), where are the roots ofξ = 0. As suggested by Fig. 1, this does not occur, since Q 0 (σ ) > 0 (except near σ = 0, discussed below).
A comment is necessary concerning the behaviour at ξ = 0. With the precise choice in (1.32), it is clear from (1.34) that for any finite value of a, Q 0 < 0 for sufficiently small σ , and alsoξ > 0 there. Thus for large but finite a, a shock will form near ξ = 0. This would slightly confuse matters, but in fact this issue is associated with the consequence at finite a that the probability density p > 0 at ξ = 0. In reality, a better initial condition would have ψ → −∞ at σ = 0, so that this region ofξ > 0 disappears, but the limit is a non-uniform one, since ξ = 0 remains a characteristic. It seems to be that the resultant shock is the cause of the necessity of treating p 0 separately.

Long-Time Evolution
The key to extending the result above is to realise that the correct way to formulate a 'continuous' distribution model is to allow p to have delta function behaviour at ξ = 0. In keeping with what the discrete model actually implies, we adopt (1.29) and thus (1.30) for ξ strictly positive, and let the probability p 0 (τ ) at ξ = 0 be finite. Thus the distribution is a Stieltjes one, and we have (1.39) where the lower limit 0+ indicates that it is in fact slightly positive (= 1 K ). Using (1.38) (which shows that ψ ≈ − 1 2 (ξ − 1) 2 near ξ = 1) together with the use of Laplace's method for the integral, we find where A is as in (1.29). The approximation in (1.31) is essentially that used in the geometric optics approximation of WKB theory (Bender and Orszag 1978), but to obtain a result equivalent to (1.18), we need the next term of the approximation. To find this, we return to (1.30), but now written in the form and we write ψ = 0 + 1 K 1 + . . .; at leading order we regain (1.31); selecting the steady solution 0 = ξ ln 1 ξ + 1 − 1, the equation for 1 is easily solved to find 1 = − 3 2 ln ξ , so that the physical optics approximation for p is This of course looks suspicious at low ξ , but it seems to be all right provided we take In fact, we now derive the equation for p 0 by takinġ is just Stirling's (rather good) approximation to 1 n! in (1.16) when n = 1. Thus the continuous probability density function approximation to p n does rather well all the way down to n = 1.
The implication of this discussion is that when a continuous model indicates a very small population being maintained for a significant time, then in practice the population will become extinct. The importance of this lies in the fact that it is not uncommon for population models which exhibit oscillations to have precisely this property, and suggests that where the oscillations are the point of focus, the continuous models are in essence incorrect.
A second conundrum which relates continuous models to low population densities is what I will call the frogspawn problem. Frogs produce thousands of eggs (e. g., Beattie 1987), but only a few of these survive to become adult frogs (Calef 1973). If we want to write a continuous model for a frog population which describes the production of thousands of eggs and their maturation as tadpoles to adult frogs, we need to find a way in which a small number of the original eggs can survive. There need to be very few, but importantly these few must not go extinct; how can that be? In a later section, I will describe one possible answer to this query.

Atto-Foxes
An enduring issue in population dynamics is what has been called 'the atto-fox problem' (Lobry and Sari 2015). The origin of this term lies in a model suggested by Anderson et al. (1981) to account for the fact that rabies outbreaks tend to recur, and they explained this by showing that oscillations can occur in the model. The work was extended by Murray et al. (1986) to account both for the recurrence and also the spread of the disease by including a term allowing for the 'diffusion' of wandering rabid foxes. An account of the model is given in his book (Murray 1989, chapter 20), since re-published in second and third editions. The epidemiology of fox rabies is described by Bacon (1985) and Toma and Andral (1977) for example. The virus is expressed in the saliva and transmitted by biting. Models of rabies dynamics continue to attract attention (e. g., Liu et al. 2017). However, the simple early models of Anderson and May and of Murray suffer from a defect, which is that the minima of the infected population reach levels which are so low as to imply extinction of the virus. Mollison (1991, p. 31) severely criticised ('this is incredible') the continuous model on two counts, the second of which is the inability of a continuous model to allow populations to become extinct, despite reaching values (in the rabies case) of one atto-fox (10 −18 of a fox) per square kilometre.
Mollison's advice was that it is essential to use a stochastic model instead of a continuous one, with the consequence of extinction. One way in which local extinction can be circumvented is through spatial heterogeneity. Most simply, if a population is distributed between relatively isolated regions with some contact, then its eradication in one region can be overcome by leakage from a neighbouring patch. For example, measles in the UK in pre-vaccination times (before 1965) was endemic in London but was transmitted to smaller communities (below the critical size necessary to maintain the endemic state) by spatial transmission (Grenfell et al. 2001;Korevaar et al. 2020), in much the same way as contraction of the heart muscle is enabled by spread of the pacemaker activity of the sino-atrial node cells to the excitable cardiac tissue cells.
The simplest way to think about spatial transmission is by the addition of a diffusion term. Of course, this only describes nearest neighbour contacts and is not suitable to describe local outbreaks caused by long-range transmission of 'sparks' (Grenfell et al. 2001), for which a spatial convolution integral would be more appropriate, but it will serve for the present discussion. It is well known that the addition of diffusion to an oscillatory system leads to periodic travelling waves. However, if the minima of the oscillations are so low as to promote extinction, then, as for measles or the heart, a local outbreak will cause the propagation of a solitary travelling wave Here I want to explore a different possibility, which is that there is something fundamentally lacking in such models, and that is the concept of a reservoir. Let me illustrate this by reference to a simple (dimensionless) population model where ε 1 represents the idea that the timescales for growth and decay of the population f are much smaller than that of the slowly varying nutrient source s. We shall see examples built around (2.1) in the following section. If s < 1, then the population plummets and will become extinct in practice. A reservoir allows for a second state g and the transfer f → g is enabled at low values of s; basically, the population hibernates until s increases above one, at which point the awakening g → f is enabled. This is somewhat similar to the concept of a refuge (Sih 1987), but it is not quite the same.
In the case of rabies, I suggested (Fowler 2000) that a resolution of the issue could be found in the distribution of infection times. SIR-type models (with a rate of recovery of the infected population I of r I ) are equivalent to assuming an exponential distribution Fig. 3 Solution of the delayed logistic equation (3.1) at α = 2.5. The asymptotic limit of large α is indicated by the flat minimum phase (the minimum of x is approximately 0.00158) of recovery times R(a) = 1 − e −ra , where a is the 'age' of the infection, so that R (a) = re −ra is the recovery time probability density (Fowler and Hollingsworth 2015). But such an assumption is not commonly realistic; if one instead assumes a fixed disease period, then the SIR model reduces to a differential-delay equation (Soper 1929). A more general assumption is that of a gamma distribution (Fowler and Hollingsworth 2015) of the form which provides a gradual change from the exponential to fixed period distributions as γ increases from 1 to ∞. In fact, experimental work in the 1960s (Parker and Wilsnack 1966;Steck and Wandeler 1980) indicated that though most rabid foxes would have incubation periods of about a month, around 10% of cases would survive for up to six months; in effect these more resistant foxes act as a reservoir for the virus. The point is that it is the small ratio of incubation time to population growth time which causes the exponentially small minima to occur.

Boom-and-Bust Dynamics
Oscillations in continuous population models which have extremely low minima occur in many other situations. One simple but remarkable example is the delayed logistic equation (Hutchinson 1948), which can be written in the forṁ In this dimensionless form, α is the ratio of the delay to the specific growth rate, and large values cause periodic solutions to occur with extremely small minima. What is remarkable is that oscillations only occur for α > 1.57, but already for α = 2.5 the asymptotic limit is visible. This is shown in Fig. 3. The minimum of the oscillations is approximately x min ≈ α exp(−e α + 2α − 1) (3.2) (Fowler 1982) and is already of order 10 −3 when α = 2.5. While the model itself is now largely only of academic interest, it is closely related to a simple model of the immune response to an infecting antigen presented by Dibrov et al. (1977a, b), which was analysed by Fowler (1981). In its simplest dimensionless form, the model is given bẏ where a and g denote antibody and antigen densities, respectively. The interpretation of the model is fairly straightforward: the infecting antigen grows but is removed by the antibodies, produced here through stimulation of the humoral immune response, the maturation time of which produces the delay in the response. If then the antigen grows unboundedly, otherwise oscillations occur, and these are severe if the delay is large. Specifically, if α and β are large (and thus necessarily κ is small), then we have a ≈ g 1 , and we regain the logistic delay equation. The behaviour of the variables is similar to that shown in Fig. 3. The immune response time is of the order of days, and much larger than a common infective (e. g., viral) growth time, and it is because of this in the model that the minima in the oscillations are attometric in scale. The immune system is a good deal more complicated than indicated in (3.3), but can be modelled in similar fashion, usually with a continuous model (Perelson 2002), and sometimes with delays (Lee et al. 2009;Rundell et al. 1998) or without (Yan et al. 2016). Commonly extinction occurs, in the sense that viral populations become very low, with the assumption that stochastic elimination occurs at these levels (Yan et al. 2016), but this is not always the case: HIV is one viral disease where an endemic state is maintained for a long period (Perelson 2002), and the herpes virus establishes itself in the body by entering a dormant state which allows for further outbreaks (Nicoll et al. 2012).

Microbial Growth
Microbial growth models are another source of oscillatory dynamics. A particularly simple example is a model due to Omta et al. (2013), who were interested in oscillations in oceanic calcifiers as a possible cause of periodic ice ages. Their model can be written asḂ  Very similar types of model have been used in describing oscillations in glycolysis (Goldbeter 1996), and plankton blooms (Huppert et al. 2005;Mahadevan et al. 2012). I represents a rate of input of nutrient to the system. It can be noted that if I = 0, (3.5) is just an SIR model, in which nutrient represents susceptibles, and the biomass masquerades as the infected. In the present case, the non-zero supply allows recovery, much as the increasing fox population allows oscillation in the rabies model. By defining the non-dimensional variables and parameter where the small parameter ε is a measure of the ratio of the nutrient supply rate to the biomass growth rate, and the nature of the model is easily understood by writing b = e θ , whence we haveθ where the potential V is defined by When ε 1, this is a slowly decaying nonlinear oscillator, and if the 'energy' E = 1 2θ 2 + V (θ ) is large, then the biomass b exhibits typically spiky 'boom-and-bust' oscillations. The form of these is shown in Fig. 4. Fowler (2014) extended this model to describe competing microbial populations (heterotrophs and fermenters). Denoting the heterotroph and fermenter populations as h and f , his model is, in dimensionless form, (3.10) here s and c are two different forms of organic carbon. The heterotrophs can utilise both, whereas the fermenters can only use the s-form, but produce the c-form, which is preferentially used by the heterotrophs. The system thus has the form of an activator-inhibitor system, and, as shown in Fig. 5, boom-and-bust oscillations occur in conditions of starvation i. e., for small ε. For ε = 0.2, the minimum of h is ≈ 1.9 × 10 −4 , but for ε = 0.1, as shown in the figure, it is ≈ 1.8 × 10 −15 . When it is reduced slightly further to the practically estimated value of ε = 0.07 , the microbial dimensionless minima are of the order of 10 −32 , and extinction beckons. Given a common bacterial loading of 10 9 cells g −1 (Kirchman 2012, p. 9), this last value corresponds to levels of 1 yocto-cell per gram, thus transcending Mollison's atto-fox per km 2 . In all these examples, extinction looms, and the oscillations are suspect. I want to suggest that extinction can be avoided in reality by the presence of a reservoir. This concept resembles but is distinct from that of a refuge for prey in predator-prey models (e. g., Sih 1987;Haque et al. 2014;Balaban-Feld et al. 2019), although the introduction of that concept was aimed at providing stability for the otherwise structurally unstable Lotka-Volterra model. The purpose and function of a reservoir in the present case is quite different. In the case of the rabies virus, a reservoir could take the form of an endemic host, which could even be the long-lived foxes themselves. But, as pointed out by Sterner and Smith (2006), there are actually many different reservoirs amongst mammals for the rabies virus and its variants. From the point of view of a continuous model, all that is necessary is that there is at least one reservoir where viral extinction does not occur.
In the case of microbes, the existence of bacteria in a dormant state is well recognised (Kaprelyants et al. 1993;Lennon and Jones 2011;Hoehler and Jørgensen 2013). Fowler and Winstanley (2018) suggested a simple model to describe the switch to dormancy; their model can be written in dimensionless form as (3.11) which generalises (3.7). Here a represents the dormant state, and p and q are switching functions which switch on (q) and off ( p) at low nutrient (c) or active biomass (b) levels. This model allows self-sustained oscillations when ε is small, again of boom-andbust type, and although generally they have less extreme minima than the yocto-cell example, effective extinction of the active biomass can occur; the difference is that the dormant bacterial population remains viable when that happens, providing a nursery for bacterial growth when the environmental stress is reduced. The reduction of b to tiny levels does not matter. This suggestion of a latent reservoir is one that can come to the rescue of continuous models when they indicate extinction. One circumstance where a reservoir may not exist is in a viral infection. Many viral infections may be completely eliminated (Rundell et al. 1998;Yan et al. 2016), although there are some where an endemic or latent state is established (Perelson 2002;Nicoll et al. 2012).

Frogspawn
Finally I want to consider another problem of small numbers and whether it can be catered for in a continuous model. Many species of plants and animals reproduce by means of the production of thousands, or even millions, of eggs or seeds. Fish provide one example (e. g., Pope et al. 2010), and trees another (Greene and Johnson 1994); helminths, discussed in the conclusions, provide another. In such cases, very few of the offspring survive to adulthood, and the question is: why? I will refer to this as the frogspawn problem, as frogs provide another example of such extreme fecundity. The population biology of frogs (it should be noted that there are many different species) has been studied by many authors (Berven 1990;Friedl and Klump 1997;Heyer et al. 1975;Smith 1983;Travis et al. 1985). Frogs produce thousands of eggs (Gibbons and McCarthy 1986), some of which later hatch to tadpoles, and still fewer of these make it to adulthood. Roughly speaking, a single adult frog in a stable population ought to produce a single offspring. How can this be?
The reason I find this perplexing is that the normal predation rate of a population F would be ∝ F (assuming plenty of predators) so that if the population becomes very low, we arrive at the previous conundrum: why does extinction not occur?
There is in fact a simple possible answer to this problem. Let us consider a population of adult frogs F, and suppose that F n is the frog population measured at intervals of T = 1 year, thus F n = F(n T ). If each female produces N eggs per year (N ∼ 1000 year −1 ), and the adult death rate is μ (year −1 ), then the year on year change in the population would be where M n (frog year −1 ) is the total number of eggs which survive predation and metamorphose to adults. In the absence of predation, M n = 1 2 N F n . The death rate is enhanced for eggs and tadpoles by predation (Waller 1973)); the reason so few tadpoles survive to adulthood is that they are consumed (even by themselves). So the principal loss term is not natural death but juvenile predation. And here is the idea: when the population levels become low, the predators do not find the prey so easily. We can for example model egg and tadpole predation during the year as a term −k M 2 , where M = 1 2 N F n initially, and There is in fact experimental evidence that this description may be reasonably accurate (see Brockelman 1969, figure 4). A suggestive continuous version of (4.1) iṡ where the overdot denotes differentiation with respect to time, and (4.4) A typical value of μ is 1 year −1 (Berven 1990). With N ∼ 10 3 year −1 , evidently δ 1 (this is the frogspawn problem), but infant predation leads to a stable population F ≈ F 0 , whose size is actually nothing to do with the egg production rate.
This of course is hardly a new idea and is the basis for Holling's type 3 predator response to prey density (Holling 1959(Holling , 1961(Holling , 1973, in which the individual predator's consumption rate as a function of prey density is S-shaped, and (for example) quadratic at low prey densities. This same idea was used in the spruce budworm model of Ludwig et al. (1978), where the birds' predation of budworm larvae is described thus by Murray (1989, p. 5): 'For small population densities …, the birds tend to seek food elsewhere…'. The motivation leading to (4.3) has the same effect as the saturating term in the logistic or Verhulst equation: whose right-hand side has the same unimodal shape as (4.3).
It should be noted that Verhulst's (1845) suggestion of the logistic equation 1 was not based on any process description, but on a wish to describe the weakening (affaiblement) of reproduction in the presence of limited resources. The present suggestion of a saturational model is based on quite different considerations.
The resolution of the frogspawn conundrum is simply due to the nonlinear egg predation rate, which allows an equilibrium to occur no matter how many eggs are produced, and whose size depends on the predation coefficient k. The surprising thing is that one normally teaches the Verhulst equation as a response to limited resources: the specific fecundity rate decreases with population size, but here the same effect is due to a quite different mechanism.

An Age-Structured Model
Of course, (4.3), while suggestive, is rather crude. A more subtle approach is to consider an age-structured model in which the amphibian density f (t, a) depends on both time t and age a. For small a, f represents eggs, then tadpoles, and for t > T , say, adult frogs. The total amphibian density is then (4.6) The units of f are taken to be Am psf −1 year −1 , where Am means amphibians, y is years, and psf is pond square foot. This last unit of area is by analogy with the Jones site model for spruce budworm outbreaks, where the larval density was measured as individuals tsf −1 , where tsf means ten square foot of susceptible branch surface area (Jones 1979;Hassell et al. 1999).
We pose the following model for f , in which the subscripts denote partial derivatives: f t + f a = −r (a) f 2 , (4.7) by analogy with (4.3). The boundary and initial conditions are taken to be where the renewal equation for f 0 is taken to be (4.9) indicating that mature adult frogs lay N eggs per year. The time T is the age of sexual maturity, commonly 1-2 years (e. g., Friedl and Klump 1997;Berven 1990). The units of r are Am −1 psf, and N has units Am Am −1 y −1 = y −1 , or eggs per frog per year.
If we define then the solution is and for t > T the renewal equation gives the integral delay equation for f 0 : . (4.12) Let us consider the form of R(a). The predation rate r is a rapidly decreasing function of a for a < T , so that R monotonically increases to a plateau at R(T ) =R, say. Thereafter R will increase slowly, and if we suppose all frogs die of senescence by age A, say, then R → ∞ as a → A. In this case the second term in the square bracket is zero for t > A.
A natural scale for R is thusR, and it is now convenient to scale the variables as (4.13) and this leads to the equation for φ, , t > α, (4.14) where = N T 1, (4.15) and the age structure is given by . (4.16) Because of the large value of , it is relatively easy to solve (4.14). We begin with an example, and consider first the steady state. R is an increasing function, with R(1) = 1 and R → ∞ as t → α. As an illustration, suppose (4.17) Fig. 6 The dimensionless age distribution given by (4.21), with = 2000 and R(a) = α(α − 1)a α(α − 1) − a(a − 1) , α = 3 (corresponding to N = 1000 y −1 , T = 2 y, A = 6 y). These are typical values for Irish frogs (Gibbons and McCarthy 1984). For visibility the range of is not shown, but at a = 0, ≈ 1530 the steady solution of (4.14) is then given uniquely by and a uniform approximation for the age distribution is (4.21) which shows that the distribution descends sharply from O( ) when a 1 to O(1) when a ∼ 1. An illustration of the resulting dimensionless age distribution is shown in Fig. 6.

Conclusions
As regards atto-foxes and yocto-cells, we suggest that the resolution of this longstanding issue may be that in practice, vanishing populations seek refuge in a safe haven, whether it be in dormancy or as an endemic remnant in another host reservoir. One example is the ability of bacteria to remain viable in the most inhospitable places (for example deep in the Earth) for an extremely long time (Hoehler and Jørgensen 2013). Plant seeds are another example: metabolic processes essentially shut down until they are stimulated to re-emerge (Bewley 1997;FitzGerald and Keener 2021).
The suggested resolution of the frogspawn problem, whether it be the logistic-type equation (4.3) or the mildly more interesting age-dependent model, seems straightforward, but it raises another issue. If the reduction of the large numbers of eggs is due to an effectively quadratic predation rate, what then is the point of the large number? We can see from (4.21) that it does not really matter how large N and thus is. It is possible that this is controlled by actual space limitation, but also, production of a small number of eggs presumably requires parental care, which is not so easily available for the predator-prone frog, and the large number simply indicates this (Davis and Roberts 2005).
There is a related issue in continuous modelling which arises in a classical model for human infection of the roundworm Ascaris lumbricoides. Roundworm infection is endemic in low-income populations with poor sanitation in tropical countries and is one of a number of 'neglected tropical diseases' which are a focus of much international interest (Holland 2013;Hollingsworth et al. 2015). The classic model to describe the infection was put forward by May (1985, 1991) and takes the essential formL = r M − μ 2 L, M = ν 0 L − μ 1 M. (5.1) Here M is the adult worm burden in the human small intestine, typically of the order of 10-20 per human, while L is the number of mature eggs in the environment. The lifetimes of eggs and adults are respectively taken to be μ −1 2 ∼ 28-84 days, and μ −1 1 ∼ 1-2 years. Actually there is an issue here already, because viable Ascaris eggs in latrines have been found with ages of the order of up to 15 years (WIN-SA 2011, Chris Buckley, private communication). Leaving that aside, (5.1) is of course linear, but nonlinearity is introduced by a unimodal dependence of the recruitment rate ν 0 on M: at low M this is due to the increasing probability of having male and female worms present in the human, and at high M because of the reduction in fecundity due to crowding. This nonlinearity allows for a stable endemic population.
Aside from the issue of the egg lifetime, there is an issue concerning the recruitment rate ν 0 . Anderson and May (1991) effectively avoided estimating this by using measured values of the basic reproduction rate R 0 = ν 0 r μ 1 μ 2 (5.2) in the range 1-5. Adult worms produce up to 2 × 10 5 eggs per day, so that in a village community of 100 people, we might have r ∼ 10 7 d −1 . It then turns out that the transmission coefficient (i. e., rate of uptake of mature eggs in the environment) is ∼ 10 −10 d −1 , something less than a nano-egg per human per day (Fowler and Hollingsworth 2016). We are back with the atto-fox problem, and the problem is worsened if we select a lower value of μ 2 . Actually this is more of a frogspawn-type problem: huge numbers of eggs only result in a small number of adults. While the latter can be understood by crowding effects in the small intestine, it does not explain why the basic reproduction rate is so low. We do not offer a resolution of this conundrum here, but the frogspawn discussion suggests that a more detailed examination of egg survival (which is assumed as a linear decay rate in deriving (5.1)) might be worthwhile.