Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Bansal, Neha
Master of Science - MSc
University of British Columbia
Reproductive value is the relative expected number of offspring produced by an individual in its remaining lifetime. For time-homogeneous population models, there are two existing meth- ods for obtaining the reproductive value function: Perron-Frobenius (PF) method and renewal method. The PF method proceeds by looking for an invariant function or an eigenvector and the renewal method proceeds by fixing a type (state of an individual) and counting the number of offspring collected there. Both methods fail for time inhomogeneous population models (those with time dependent parameters), therefore other methods are needed to compute the reproductive value in that case. Towards that goal, we present two new characterizations of this function: for a general class of critical age-structured population models, we show that the asymptotic ratio of the survival probabilities, from different initial states, is equal to the initial ratio of reproductive values, and that in relation to a specific coupling of the size-biased model, both are equal to the limit of the inverse ratio of total reproductive values. To provide context and background, we also discuss two defining properties of reproductive value, namely time-invariance in deterministic models and the martingale property in stochastic models, for some commonly used models. These findings contribute to a deeper understanding of reproductive value functions and their implications in critical population models.
Reproductive Value for Branching Population ModelsbyNeha BansalB.Tech., Indian Institute of Technology Roorkee, 2015A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)August 2023© Neha Bansal, 2023The following individuals certify that they have read, and recommend to the College ofGraduate Studies for acceptance, the thesis titled:Reproductive Value for Branching Population Modelssubmitted by Neha Bansal in partial fulfilment of the requirements of the degree ofMaster of Science in MathematicsExamining Committee:Dr. Eric Foxall, Assistant Professor, I. K. Barber Faculty of SciencesSupervisorDr. Rebecca Tyson, Professor, I. K. Barber Faculty of SciencesSupervisory Committee MemberDr. Michael Noonan, Assistant Professor, I. K. Barber Faculty of SciencesSupervisory Committee MemberDr. Loic Markley, Associate Professor, School of EngineeringUniversity ExamineriiAbstractReproductive value is the relative expected number of offspring produced by an individual inits remaining lifetime. For time-homogeneous population models, there are two existing meth-ods for obtaining the reproductive value function: Perron-Frobenius (PF) method and renewalmethod. The PF method proceeds by looking for an invariant function or an eigenvector andthe renewal method proceeds by fixing a type (state of an individual) and counting the num-ber of offspring collected there. Both methods fail for time inhomogeneous population models(those with time dependent parameters), therefore other methods are needed to compute thereproductive value in that case. Towards that goal, we present two new characterizations ofthis function: for a general class of critical age-structured population models, we show that theasymptotic ratio of the survival probabilities, from different initial states, is equal to the initialratio of reproductive values, and that in relation to a specific coupling of the size-biased model,both are equal to the limit of the inverse ratio of total reproductive values.To provide context and background, we also discuss two defining properties of reproductivevalue, namely time-invariance in deterministic models and the martingale property in stochasticmodels, for some commonly used models. These findings contribute to a deeper understandingof reproductive value functions and their implications in critical population models.iiiLay SummaryThe notion of reproductive value provides a measure of an individual’s and a population’s po-tential for reproduction at any given time. In this study, we explain the techniques used tocalculate the reproductive value in various population models, which are either predictable (de-terministic models) or involve some randomness (stochastic models). One significant findingis that when the rates of birth and death in a population stay the same over time and balanceeach other out, the reproductive value is constant for deterministic models and is constant onaverage for stochastic models.In population models where each individual has an average of one offspring, we present anovel finding: as time t becomes large, the ratio, from different initial states, of the probabilityof surviving until time t tends to the ratio of the reproductive values of those initial states.Especially, initial states with higher reproductive value have relatively longer-lived lineages.ivContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Chapter 2: Reproductive value for deterministic population models . . . . . . . . 72.1 What is reproductive value? . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1.1 Perron-Frobenius (PF) theory . . . . . . . . . . . . . . . . . . . . . . 82.1.2 Renewal Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Linear Population models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.1 Age Structured Model . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.2 Multi Type Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.3 Non-linear Population Models . . . . . . . . . . . . . . . . . . . . . . . . . . 232.3.1 Periodicity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.3.2 Uniform Non-linearity . . . . . . . . . . . . . . . . . . . . . . . . . . 262.3.3 State Dependent Non-linearity . . . . . . . . . . . . . . . . . . . . . . 29Chapter 3: Stochastic Population Models . . . . . . . . . . . . . . . . . . . . . . . 303.1 Branching Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.1.1 Discrete Time Scenario . . . . . . . . . . . . . . . . . . . . . . . . . . 323.1.2 Continuous Time Scenario . . . . . . . . . . . . . . . . . . . . . . . . 35v3.2 Multi-type Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.3 Age structured population model . . . . . . . . . . . . . . . . . . . . . . . . . 42Chapter 4: New characterizations of reproductive Value . . . . . . . . . . . . . . 464.1 Reproductive value and survival probability . . . . . . . . . . . . . . . . . . . 484.2 Size Biasing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.2.1 Martingale change of measure . . . . . . . . . . . . . . . . . . . . . . 564.2.2 Spine process construction . . . . . . . . . . . . . . . . . . . . . . . . 594.2.3 Age Structured Population Model . . . . . . . . . . . . . . . . . . . . 614.3 Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.3.1 Discrete time random walk . . . . . . . . . . . . . . . . . . . . . . . . 674.3.2 Continuous time size biased binary branching process . . . . . . . . . 684.4 Reproductive value and martingale . . . . . . . . . . . . . . . . . . . . . . . . 70Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87viList of FiguresFigure 2.1 Invariant function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 2.2 Invariant measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Figure 2.3 Age structured population model . . . . . . . . . . . . . . . . . . . . 18Figure 3.1 Branching tree diagram. . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 4.1 Coupling of two size-biased critical branching processes . . . . . . . . 74viiAcknowledgementsI would like to express my gratitude to Dr. Eric Foxall for his consistent guidance and encour-agement during my graduate studies. His unwavering support and comprehensive feedbackmade my time as his student truly enjoyable. I am extremely thankful for his remarkable atten-tion to detail, which greatly facilitated the challenging task of writing. Additionally, I wouldlike to extend my appreciation to Dr. Rebecca Tyson and Dr. Michael Noonan, members ofmy supervisory committee, for dedicating their time to review my thesis and provide valuablefeedback.Furthermore, I would like to acknowledge the Tyson Lab (Tylab) research group for theirsupport and valuable insights throughoutmy tenure as a graduate student at UBCO. I am gratefulto my roommate, Shambhavi Singh, for keepingme sane. The laughter andmemories we sharedwill forever hold a special place in my heart and memory.viiiChapter 1IntroductionPopulation biology is the study of population dynamics, such as the growth and decline of pop-ulations, interactions among species, long-term behavior, and the effect of external factors. Tomodel biological and ecological systems, various mathematical models are used. These modelsare often structured models, i.e., individuals have a type (state) that changes over time and af-fects their behaviour. Individuals with different types have varying abilities to contribute to thefuture success of the population. As long as the parameters that define individual behavior areconstant over time (which is nevertheless a strong assumption), we can usually define, in at leasttwo different ways, a reproductive value function that quantifies these relative contributions, asa function of type. It would be interesting to have these functions for a larger class of models,in particular those with time-dependent parameters. However, existing methods break down inthis case. Consequently, in this work, our goal is to find new characterizations of the reproduc-tive value function, with the intention that one or more of these characterizations can eventuallybe used to generalize the notion to more complicated models. The mathematical models that weconsider are broadly categorized into two classes: deterministic models and stochastic models.Deterministic population models are formulated using differential equations for studyingcontinuous-time behavior and using difference equations for studying discrete-time behavior.If the system’s differential equation can be solved, it is possible to predict the exact state ofthe system; otherwise, approximations or steady-state solutions are used to describe long-termpopulation dynamics. More details in the mathematical formulation of deterministic models are1added by accounting for interaction among individuals, which are known as non-linear deter-ministic models. Generally, deterministic population models capture the general form of linearor non-linear relationships for an ecological system in a deterministic manner.Ecological systems are formulated as stochastic population models to include the random-ness that arises in real-world systems due to many factors, such as environmental and demo-graphic factors. This randomness can arise in the timing of individual events (demographic) orin the parameters that define the model (environmental). Some stochastic population modelscan be formulated using stochastic differential equations.The stochastic models in which individuals are dependent on each other are known as in-teracting particle processes. An example is the voter model [Liggett et al., 1999], where theopinion of an individual changes by interacting with other individuals. In this thesis, we areconsidering that individuals are independent of each other, i.e., the behavior of an individualis not impacted by other individuals. These stochastic models are known as branching particlemodels, where an individual’s birth state defines its lifetime and reproduction behavior.The branching population models in which individuals in a population follow an identicaldistribution for the lifetime and reproduction process are known as Crump-Mode-Jagers (CMJ)processes or general branching processes [Jagers, 1989]. A CMJ branching process is a stochas-tic process in which an individual is characterized by two random objects: the (non-negativereal) lifespan and the offspring measure, which is a counting measure on R+. All individualsare independent and follow these two identical distributions, although the reproduction processand lifetime can be dependent on each other for an individual.The deterministic version of a stochastic population model can be derived by setting the rateof change in the equation equal to the ‘drift’, i.e., an instantaneous rate of change of expectedvalue at each state. In this way, the deterministic model characterizes the mean behavior ofthe stochastic system. If we go in reverse, that is, from deterministic to stochastic, we haveto assign distributions to the timing of events such as birth and death; the simplest way to dothis is using Poisson processes with appropriate intensities such that the mean behavior of thestochastic model coincides with the deterministic model.Our main objective is to characterize, in various ways, reproductive value for various2stochastic and deterministic population models. We are focusing on obtaining the formula-tion of reproductive value for time-homogeneous population models (birth and death rates arenot dependent on time). For obtaining these characterizations, we need a martingale, i.e., astochastic process for which expected value at any time t conditioned on information at anytime s < t is equal to its value at time s, and we need a reproductive value function to get themartingale.Martingales are a key tool in probability theory for studying stochastic processes. In thisthesis, we make use of so-called additive martingales, which are obtained by integrating somefunction against the population measure and possibly multiplying by some function of time tocompensate. Martingales are (typically real-valued) stochastic processes that are ‘constant onaverage’ in a sense made precise in [Durrett, 2019]; they generalize the concept of a conservedquantity to stochastic processes. They have nice convergence properties under mild assump-tions, for example, if a martingale Xt has bounded L1 norm, i.e., suptE[|Xt|] < ∞ then ast→∞, Xt converges almost surely to a limit X with finite expectation, a precise definition isin [Durrett, 2019]. The convergence property is useful for characterizing the asymptotic behav-ior of stochastic processes. Also, good estimates are available for the extreme values of theirtrajectories such as Doob’s inequality, i.e., we can bound the supremum by the distribution ofthe process at the initial or final time. Martingales are used for adapted change of measure,in other words, if we want to understand the same stochastic process under different condi-tions or by biasing, then the quantity that translates between the original and changed measure(formally, the Radon-Nikodym derivative of the changed measure with respect to the originalmeasure, conditioned on information up to time t, as a function of t) is a martingale.Martingales are used for analyzing the long-term behavior of branching processes by es-tablishing the conditions for their convergence. In [Hong et al., 2022], authors proved the con-vergence of martingale for critical discrete time branching process in a random environment,i.e., offspring distribution is varying across generations. Martingales are widely used to findconditions for extinction and the explosion of a branching process, authors in [Li et al., 2019]established these conditions for a continuous state branching process (e.g. age structure pop-ulation models), where the birth rate is dependent on the population density and competition3among individuals are taken in account. For numerical simulation of stochastic processes, themartingale formulation is utilized for estimating approximate equilibrium behavior, i.e., sta-tionary distribution by [Belomestny et al., 2022], authors used their approach for simulatingwell-known Lotka-Volterra predator-prey system.In [Kyprianou and Palmowski, 2005], authors analyzed and proved the limiting behaviorresults for the Levy processes (Brownian motion) using martingales and Markov properties.In the study [Kyprianou, 2000] of spatial branching random walk, which is similar to spatialpopulation model in discrete time, authors used a size-biased branching process to establish theconvergence of martingales. We are using the size-biased method to compare two branchingprocesses via martingale convergence. In [Biggins and Kyprianou, 2004], martingale changeof measure is used to establish the conditions for the asymptotic behavior of general branchingprocesses (CMJ branching process). In our work, we are using a martingale to define the mea-sure for the size-biased branching process to establish the relation between martingale and thereproductive value of the population process. On a similar line, [Kolesko and Sava-Huss, 2023]generalizes the results related to the long-term behavior of multi-type branching processes usingmartingales. In [Jagers and Klebaner, 2000], reproductive value and the long-term behavior ofthe general branching processes with birth rate dependent on population density and age areexplored.In the context of branching processes, the reproductive value functions are also knownas eigenfunctions or invariant functions because of the invariant (martingale) propertyof the reproductive value. The eigenfunctions are used in [Pavliotis and Zanoni, 2022][Ojeda and Pardo, 2020] for an interacting particle system, i.e., a branching processwith interacting individuals, to analyze the long-term behavior of the system. In[Bianca and Dogbe, 2017], authors prove the uniqueness and existence of the invariant functionfor a multi-type branching process. In [Molchanov and Yarovaya, 2012], the invariant func-tion is used for bifurcation analysis for a branching random walk. In [Moharrami et al., 2020],authors derive an eigenvalue and an eigenfunction for a multi-type branching process with un-accountably many types and use it to study the probability of extinction. For a continuousstate branching process, in [Foucart and Vidmar, 2022], authors are using the eigenfunction4for analyzing a general continuous-state branching process. In [Cox et al., 2022], the eigen-function defines the large population limit of a binary branching process with interacting par-ticles, where binary branching means two or zero offspring are produced at each birth event.In [Cloez et al., 2021], long-term behavior is analyzed for a size-structured branching processusing the eigenfunction, which is similar to the age-structured branching process that we arestudying in this thesis.In Chapter 2, linear and non-linear deterministic population models are discussed in contin-uous time and discrete time. Specifically, the age-structured population model and multi-typepopulationmodels are described in detail. Here, we also discuss interactions among individuals,which arise from the dependency of birth and death rates on the population density, an individ-ual’s state, or time. Some examples of non-linear deterministic models are the predator-preymodel and infectious disease transmission models. In Chapter 3, single-type and multi-typebranching processes are discussed in discrete and continuous time, which are categorized ac-cording to the value of the Malthusian parameter [Fisher, 1999]. The Malthusian parameter isthe growth rate, assuming that the population grows exponentially. In Chapter 3, we discussa special case of the CMJ process, which is an age-structured population model in continuoustime in detail, in which birth rate and death rates are dependent on age.Chapters 2 and 3 give an introduction to established methods for estimating reproductivevalue functions and the invariance property of the reproductive value. Specifically, in Chapter 2we discuss the invariant functions for various well-known deterministic population models andin Chapter 3 we discuss the martingale properties for well-known stochastic population models.In Chapter 4 we focus our attention on the critical branching process, in which case, on av-erage, a particle exactly replaces itself in its lifetime. We establish two new characterizations ofthe reproductive value function: the first one is proof of the relation between the ratio of repro-ductive value and the ratio of survival probability, and the second one is the relation betweenthe ratio of reproductive value and the inverse ratio of martingales, for the time-homogeneousage-structured stochastic population model. The motivation for this characterization is to pro-vide an approach for defining reproductive value that could generalize to time in-homogeneouspopulationmodels, in which case birth and death rates may depend explicitly on time. That said,5the problem of finding a reproductive value function for a wide range of time-inhomogeneousmodels is still open.6Chapter 2Reproductive value for deterministicpopulation models2.1 What is reproductive value?In population genetics, reproductive value is defined by Ronald Fisher in [Fisher, 1999] as thenumber of female offspring born to a female in its remaining lifetime relative to the population.In other words, reproductive value, which is a function of age a, is equal to the average contri-bution of an individual of age a to the ancestry of future generations [Lion, 2017]. There arebroadly two methods to estimate the reproductive value of a population model. One is usingrenewal/regeneration theory and the other is using Perron-Frobenius (PF) theory. Below, wediscuss both approaches. The PF approach is essentially a linear algebra/functional analysisapproach where you look for the spectrum to have a certain property, such as a distinct largesteigenvalue that is positive along with a positive eigenvector, as well as a spectral gap. On theother hand, the renewal approach is similar to using a return map in a dynamical system, orusing the strong Markov property in a Markov chain: shrinking your observations to a subsetof the state space and keeping track only of its returns to that set.Definition 2.1.1. A Markov chain is a stochastic process in which the probability to take thenext move is dependent on the current state of the process and independent of the history of theprocess, which is also known as the Markov property. More precisely, we sayXn is a discrete-7time Markov chain with p(i, j) as probability of transition from state i to state j, if for anyj, i, in−1, ..., i0P (Xn+1 = j|Xn = i,Xn−1 = in−1, ..., X0 = i0) = p(i, j).To define the reproductive value function, in the perspective of PF theory, is to look fora invariant function and in the perspective of renewal theory, we fix a type (age 0 in an age-structured population model) and count the number of offspring collected there. In general,the reproductive value of an individual depends on its state at the time of estimation. In timehomogeneous population models, i.e., birth and death rates don’t depend on time, reproductivevalue is a function only of the state, whereas for time in-homogeneous population models inwhich birth and death rates are dependent on time, the reproductive value, if it exists, maydepend on not only the state but also time. The inherent time dependency of birth and deathrates poses a challenge for the naive definition of reproductive value, as it fails to exhibit thecrucial property of remaining constant over time.Both the renewal and PF theory definitions of reproductive value can be applied only totime-homogeneous models, the reason for which is explained in Sections 2.1.1 and 2.1.2. So, ifwe wished to extend the notion of reproductive value to time-inhomogeneous models we wouldneed a more general definition that, for time-homogeneous models, still returns the functionobtained from the PF or renewal approach. This is what motivates our first characterizationof the reproductive value function, modulo positive multiples, as the limiting ratio of survivalprobabilities starting from different initial conditions, as discussed in Chapter 4.2.1.1 Perron-Frobenius (PF) theoryA matrix A is non-negative if all the elements of A are non-negative i.e aij ≥ 0 ∀i, j and A isa positive matrix if aij > 0 ∀i, j. A non-negative or positive matrix is regular if there existsa k ≥ 1 such that Ak > 0, i.e., all elements of matrix Ak are positive. PF theory asserts thatfor a real non-negative square matrix, there exists a unique positive eigenvalue of the largestmagnitude (also known as dominant eigenvalue) and corresponding non-negative eigenvectors8which are unique in the case of positive matrices. In matrix theory, the PF theory was provedby Oskar Perron in [Perron, 1907] for regular positive matrices and by Georg Frobenius in[Frobenius et al., 1912] for regular non-negative matrices.For real non-negative and regular matrices, Perron-Frobenius theory states that there is adominant real eigenvalue λpf ≥ 0 with multiplicity 1 and associated non-negative left andright eigenvectors denoted by v and w respectively which are unique up to positive multiples.The right eigenvector w is the reproductive value vector. A common normalization is to take1Tv = 1 and wTv = 1.In the context of population models, the matrix A is a reproduction matrix[Fujiwara and Caswell, 2001] i.e element aij of matrix A represents the number of offspringof type j produced by individuals of type i in discrete time setting. Here, type can be a de-mographic or spatial characteristic of the individuals such as age, size, life stage, health status,etc. Note that, A is always a non-negative matrix for population models but it may not alwaysbe a regular matrix. Therefore, the use of PF theory is limited to cases where the reproductionmatrix is regular.If we have a non-negative and regular reproduction matrix A then the left eigenvector vw.r.t λpf of A gives the asymptotic population distribution, and the right eigenvector w givesthe reproductive values. In particular, wi is the reproductive value of a type i individual andvi is the asymptotic fraction of type i individuals. The PF eigenvalue λpf gives the asymptoticgrowth rate i.e in long-term if λpf > 1 then population is growing, if λpf < 1 then population isdecaying and if λpf = 1 then population reaches an equilibrium distribution as t → ∞, whichis the right eigenvector of A with respect to λpf. Furthermore, when λpf = 1, total reproductivevalue is constant over time which is equal to the inner product of the eigenvector w and thepopulation vector at any time t.Theorem 2.1.1. [Li and Schneider, 2002] Let A be a non-negative regular matrix andX(n) isgiven by,X(n) = X(0)An,9where X(0) is an arbitrary, fixed initial value, then for large n,X(n) ∼ λtpfwTX(0)v, (2.1)where λpf is dominant eigenvalue of A, w is right eigenvector and v is left eigenvector w.r.t λpf.Age structured population model is a typical example, in which case reproduction matrix Ais also known as Leslie matrix [Leslie, 1945], where ai1 denotes the fecundity of individuals ofage i, for all i, ai,i+1 denotes the number of individuals that survive from age i to age i+ 1 andremaining elements are 0 as shown below.A =f1 s1 0 . . . 0 0f2 0 s2 . . . 0 0. 0 0 s3 . . . 0. 0 0 0 . . . 0. 0 0 0 . . . 0fI 0 0 0 . . . sI(2.2)More generally, PF theory can be extended for non-negative compact operators by theKrein-Rutman theorem [Krein and Rutman, 1948]. A compact operator is a linear operatorT : X → Y , whereX and Y are normed vector spaces such that T maps bounded subsets ofXto relatively compact subsets of Y . The reproduction matrix from the age-structured model isfinite rank, and thus a compact operator. Krein-Rutman theorem states that for a compact andnon-negative operator T , there exists a dominant eigenvalue with an eigenvector. This theoremis used for defining reproductive value when the reproduction matrix is infinite dimensionalsuch as when spatial movement is included in the model.2.1.2 Renewal TheoryLet τ1, τ2, ... be random independent exponential variables with parameter λ. Let Tn =∑ni=1 τifor n ≥ 1 and T0 = 0 and define N(t) = max{n : Tn ≤ t}. Here, N(t) is the number ofevents that occur in the time interval [0, t], where N(0) = 0 and τi is the time interval between10the occurrence of i − 1 and ith event. N(t) is a counting process that has Poisson distributionwith rate parameter λ. Note that, the time between event occurrences are independent of eachother but follow the same exponential distribution. In other words, the process gets renewedor it starts afresh after every event occurrence. A renewal process is a generalization of thePoisson process where time intervals are independent and identically distributed (i.i.d) randomvariables but may not have an exponential distribution.In the context of population models, renewal occurs through the cycle of birth and death.Hence, by tracking the birth events we can keep an account of the survival of the population.Now, assuming individuals are independent and follow similar distribution for reproductionprocess and survival probability, the renewal process is like bookkeeping of population changesby tracking birth events [Lotka, 1939a].To use the renewal theory definition, it is assumed that a newborn individual evolves intime independent of its parent as well as other individuals but follows the same distributionfor the reproduction process and survival probability. This assumption is the basis for defininggeneral branching processes also known as CMJ processes. We will discuss the age structuredbranching process in Chapters 3 and 4, which is a special case of CMJ processes.The existence of a renewal point in population processes makes the analysis of long-termpopulation dynamics less complicated. Here, we are assuming that the renewal point is a recur-rent point for the population process, which means that the process needs to hit the renewal/re-current point infinitely many times to survive. Using renewal points, we can partition the pro-cess by renewal points in time. The end and start point of a partition is the renewal point. Theseexcursions are i.i.d because future birth events of an individual are independent of its past butfollow the same reproduction distribution. This property of the renewal process is similar tothe strong Markov property (see Chapter 1, Theorem 1.2 in [Durrett, 1999]), which states thatfor a Markov process probability of moving to the next state depends on the current state ofthe process and it is independent of the history of the process [Durrett, 1999]. Note that, theuse of renewal theory is not limited to CMJ population processes only, it can be used for anypopulation process whose reproduction kernel has a renewal point. In this thesis, however, weare focusing on CMJ processes.11Further, if the birth state of all the individuals is the same irrespective of the parent’s statethen reproductive value is defined as the expected number of births of an individual in its remain-ing lifetime. The common birth state that signifies the occurrence of a birth event or renewal isalso known as the renewal state of the process. For example, in age structured population mod-els, age 0 is a renewal state. Note that, to survive for a long time, the population process needsto hit the renewal state many times. This resembles the concept of recurrent state in Markovprocesses [Durrett, 1999].Let’s understand the approach and ingredients required to define reproductive value by re-newal theory definition, following the treatment given in [Athreya and Ney, 1982]. Consider apopulation process in discrete time which is defined on a measurable space (Ω,F) with statespace as Ω and a sigma algebra F of subsets of Ω. Let K : Ω × F 7→ R+ be a reproductionkernel, which is analogous to the reproduction matrix from PF theory. Kn is the nth iterate ofK for n ≥ 0 and satisfiesK0(x,E) = χE(x), K1(x,E) = K(x,E) (2.3)Kn+1(x,E) =∫ΩKn(y, E)K(x, dy), n ≥ 1 (2.4)for all x, y ∈ Ω and for all E ∈ F , where χE(x) is the indicator which is 1 if x ∈ E else 0. Incontext with population models, Kn(x,E) is the number of nth generation offspring of type inset E produced by an individual of type x.Now, if there exists a state ∆ ∈ Ω such that,∑nKn(∆,∆)sn =<∞ for 0 ≤ s < R=∞ for s ≥ R(2.5)then ∆ is a recurrent state or renewal state [Athreya and Ney, 1982]. Here, 1/R acts as thedominant eigenvalue in PF theory, i.e., λpf. In the case of Markov chains, K is a stochastickernel whereK(i, j) is the probability to move from state i to state j in one step, and the abovesummation with s = 1 represents the expected number of visits to the recurrent state∆ startingfrom ∆. In the context of population models where K is a reproduction kernel, ∆ represents12Figure 2.1: For a population process, the reproductive value f(x) of an individual of type x isdefined by measuring the number of birth events that is reaching the state ∆ in any number ofsteps. In the case of an age-structured populationmodel, steps represent the age of the individualat the time of a birth event.the renewal point.If there is a recurrent state ∆ for this process then there exists an invariant function f(x),i.e., (Kf)(x) =∫F f(y)K(x, dy) = f(x), which is estimated by measuring the number ofindividuals of state ∆ produced by a type x individual in its remaining lifetime as shown infigure 2.1. For detailed derivation see [Athreya and Ney, 1982]. This invariant function isnothing but the reproductive value of an individual of type x. Note that, this form of invariantfunction is analogous to the right eigenvector form in PF theory with λpf = 1.Also, if ∆ is a recurrent state then there exists an invariant measure µ i.e µK(A) =∫ΣK(x,A)µ(dx) = µ(A), which is estimated by measuring the number of individuals in eachstate as t→∞, if the process starts from state∆ as shown in figure 2.2. For detailed derivationsee [Athreya and Ney, 1982]. The invariant measure is the asymptotic population distributionwhich is similar to the left eigenvector in PF theory for λpf = 1.Reproductive value definitions using renewal state and PF theory coincide. Both of themrepresent the same object which is the average number of offspring of an individual in its re-maining lifetime as shown in [Athreya and Ney, 1982].13Figure 2.2: For a population process starting with an individual of type ∆, asymptotic popula-tion distribution (µ) is estimated by measuring the proportion of individuals which reach a statein set E in any number of steps. In the case of an age-structured population model, it is theproprotion of individuals reaching the age range represented by set E.Reproductive value definition using renewal and PF theory is limited to time homogeneouspopulation models, in the case of time in-homogeneous models, birth rate and death rate aredependent on time which leads to different distributions for reproduction and survival aftereach birth event or renewal point. Therefore, the condition that excursions between successiverenewals be i.i.d is violated.2.2 Linear Population modelsIn this section, a well-known class of deterministic models is considered in which birth anddeath rates are only dependent on the age or type of an individual. These models are definedwith the assumption that the individuals are independent of each other, and their reproductionprocess and lifespan are not impacted by other individuals or other environmental conditions.We will discuss the dynamics of these models by defining the reproductive value using renewaland PF theory.A general form to represent these models in discrete time is using the difference equationXn+1 = XnA, (2.6)14where Xn is the population vector if individuals are differentiated by their state or type and Ais a reproduction matrix, which is non-negative and doesn’t change with time. For continuoustime and discrete state population models, the general form is given using differential equationssuch asdX(t)dt= X(t)A, (2.7)whereA is a constant reproduction matrix with non-negative off-diagonal terms andX is popu-lation vector at time t. In the following section, we will discuss a population model where state(age) is a continuous variable and the model is formulated using partial differential equations.2.2.1 Age Structured ModelA population can be structured with respect to many features on the basis of demog-raphy and physiology such as age, developmental stage [Manly, 1990], size, and gender[Al-Arydah et al., 2012]. In [Auger et al., 2008], authors study the population model with age-size-spatial structure in biology and epidemiology. An example of size structured populationmodel is developed in [Rose et al., 2005] for Platte thistle (Cirsium), to study the impact of sizeon the flowering of the plant. Age is one of the prominent features for structuring the populationbecause birth and death rates are significantly impacted by the age of an individual. The study ofage structured models was initiated by A.J. Lotka and A.C. M’kendrick in [M’kendrick, 1925][Lotka, 1939b], and [Smith et al., 1977].Consider a population model in a continuous time setting with age dependent birth rate β(a)and death rate µ(a). Assuming the maximum age of each individual is A and the birth age is0. Let u(a, t) denote the population density of individuals of age a at time t. Hence, the totalpopulation at time t is,N(t) =∫ A0u(a, t)da. (2.8)The initial state of the population at t = 0 is a given age distribution denoted by φ(a) fora ∈ [0, A]. The change in population density is formulated by following a system of equations15also known as McKendrick-Von Foerster [Murray, 2002] equation,∂au(a, t) + ∂tu(a, t) = −µ(a)u(a, t) (2.9)u(0, t) =∫ A0β(s)u(s, t)ds (2.10)u(a, 0) = φ(a) ∀a ∈ [0, A], (2.11)where ∂au(a, t) represents the aging component, ∂tu(a, t) denotes the change in density withtime. For any age a 6= 0, in a small time interval [t, t+h], the change in the density of the cohortthat has age a at time t (and will have age a+ h at time t+ h) is approximately−µ(a)u(a, t)h.The density of newborn individuals u(0, t) is integral of the total birth rate at time t over allages, where the total birth rate at time t for age a is given by β(a)u(a, t)ds.It is assumed that individuals are independent of each other, i.e., the birth rate and death ratedon’t depend on the lifetime process of other individuals. As explained in Chapter 1, this modelgives the mean behavior of a stochastic process in which the birth age of individuals is 0 andthey evolve in time independently from their parent and other individuals. This indicates thatthe population is renewed with every cohort of new births and can be studied independently.Therefore, age 0 is a renewal state for this model. Now, we can utilize the renewal definitionto find the reproductive value of an individual of age a as,h(a) =∫ Aap(s)p(a)β(s)ds, (2.12)wherep(a) = exp(−∫ a0µ(s)ds) (2.13)is the probability that an individual will survive till age a. To make this model more realistic,the birth rate can be defined as a function of age such that it is 0 up to the maturity age as wellas after the fertility period ends. In between, it can be a parabolic function or another type offunction in accordance with the population under consideration.Before classifying the age structured population model, we give a precise definition of the16Malthusian parameter.Definition 2.2.1. [Malthus and Winch, 1992] The Malthusian parameter is the parameter r inthe exponential equation of population growth N(t) = N0 exp(rt), where N0 is the initialpopulation size (at t = 0) and t is the elapsed time.The age structured population models are categorized into 3 classes in accordance withMalthusian parameter α for the process which is calculated by solving the following equation,∫ ∞0e−αtp(t)β(t)dt = 1. (2.14)Here, p(t)β(t) is the expected rate of offspring production at age t and e−αt is the discountingfactor such that integral over all ages is 1. Hence, the process is classified as follows,• if α < 0, then the process is sub-critical, as it will compensate for fewer birth events tomake the integral value 1.• if α = 0, then the process is critical and• if α > 0, then the process is super-critical, as it will compensate for more birth events tomake the integral value 1.Next, for a critical age structured population model, we present the property of reproductivevalue in the following lemma, which is well-known in the field, although we can not find aprecise reference. So, we included the proof here for your convenience.Lemma 2.2.1. LetH(t) =∫ A0h(a)u(a, t)da be the reproductive potential of the population attime t. If age structured population model is a critical process then the reproductive potentialof the population is a time-invariant function.Proof. Let’s estimate the rate of change of H(t),dH(t)dt=∫ A0du(a, t)dth(a)da (2.15)=∫ A0(−µ(a)u(a, t)− ∂au(a, t))h(a)da (2.16)17Figure 2.3: Age structured population process starting with an individual of age 0 at time t.= −∫ A0∂au(a, t)h(a)da−∫ A0µ(a)u(a, t)h(a)da (integration by parts) (2.17)= −[h(a)u(a, t)]A0 +∫ A0u(a, t)h′(a)da−∫ A0µ(a)u(a, t)h(a)da (2.18)= u(0, t)h(0)− u(A, t)h(A) +∫ A0u(a, t)h′(a)da−∫ A0µ(a)u(a, t)h(a)da (2.19)= u(0, t)h(0) +∫ A0u(a, t)h′(a)da−∫ A0µ(a)u(a, t)h(a)da, (2.20)where last equality holds because h(A) = 0. Next, calculating h′(a),h(a) =∫ Aap(s)p(a)β(s)ds (multiply both sides by p(a)) (2.21)h(a)p(a) =∫ Aap(s)β(s)ds (differentiate both sides w.r.t a) (2.22)h′(a)p(a) + p′(a)h(a) = p(a)β(a) (2.23)h′(a)p(a)− µ(a)p(a)h(a) = p(a)β(a) (2.24)h′(a)p(a) = µ(a)p(a)h(a) + p(a)β(a) (2.25)h′(a) = µ(a)h(a)− β(a). (2.26)18Hence dH(t)/dt becomes,dH(t)dt= u(0, t)h(0) +∫ A0u(a, t)(µ(a)h(a)− β(a))da−∫ A0µ(a)u(a, t)h(a)da (2.27)= u(0, t)h(0)−∫ A0u(a, t)β(a)da = u(0, t)h(0)− u(0, t) (2.28)= (h(0)− 1)u(0, t). (2.29)Thus, for H(t) to be invariant, we want dH(t)dt= 0, i.e., h(0) = 1 that means h(0) =∫ A0p(s)β(s)ds = 1. Hence, forH(t) to be invariant, the expected number of births by a particlein its lifetime should be one.Lemma 2.2.2. Let u(a, t) be the population density of an age structured population model andlet p(a) be as in (2.13). If PDE is initialized with u(a, 0) = p(a), the model is a critical process,and β(a) is bounded then u(a, t) = p(a) for all t > 0.Proof. Consider a cohort of population at time t0 of age a0. Let u¯(h) = u(a0 + h, t0 + h) andµ¯(h) = µ(a0 + h). Then (2.9) is,du¯dh+ µ¯(h)u¯ = 0. (2.30)Solving the above ODE we get,u¯(h) = u¯(0) exp(−∫ h0µ¯(s)ds). (2.31)Hence, the solution is,u(a, t) =p(a)p(a−t)u(a− t, 0) if a > t,p(a)u(0, t− a) if a < t.(2.32)Now, the initial condition u(a, 0) = p(a) implies that u(a− t, 0) = p(a− t). So, for a > t,u(a, t) = p(a)p(a−t)p(a− t) = p(a). In order to show that u(a, t) = p(a) for a < t, it is sufficientto show that u(0, t) = p(0) = 1 for all t > 0. So, let’s check that by expanding u(0, t) using19the general solution of u(a, t) from (2.32),u(0, t) =∫ A0β(a)u(a, t)da (2.33)=∫ A∧t0β(a)p(a)u(0, t− a)da+∫ AA∧tβ(a)p(a)p(a− t)u(a− t, 0)da (2.34)using the initial condition u(a− t, 0) = p(a− t) (2.35)=∫ A∧t0β(a)p(a)u(0, t− a)da+∫ AA∧tβ(a)p(a)da, (2.36)where A ∧ T = min(A, T ). The equation (2.36) has the form of the integral equation (1.1) inreference to [Feller, 1941] with f(a) = β(a)p(a) for 0 ≤ a < A and f(a) = 0 for a > A,g(a) =∫ AA∧t β(a)p(a)da and u(t) = u(0, t). In that reference, it is shown that there exists aunique solution u(t) when f and g are non-negative and bounded functions in any boundedinterval. Note that, f(a) and g(a) are non-negative as birth rate β(a) is non-negative for alla ∈ R+ and p(a) is the survival probability which is a non-negative quantity. Assuming thatβ(a) is bounded implies that f(a) and g(a) are bounded functions. So, it suffices to check thatu(0, t) = 1 for all t > 0 works. Let’s substitute u(0, t− a) = 1 in (2.36),u(0, t) =∫ A∧t0β(a)p(a)da+∫ AA∧tβ(a)p(a)da (2.37)=∫ A0β(a)p(a)da (2.38)= 1. (2.39)This implies that u(0, t) = 1 for all t > 0 is the solution of (2.33). Hence, for both a > t anda < t, u(a, t) = p(a).Example: Consider an age structured population model with birth rate β(a) = 4a and deathrate µ(a) = 2a. Here, survival probability is p(a) = exp(− ∫ a02sds) = exp(−a2). Supposethe maximum age is 50 years then the reproductive value of an individual of age a is,h(a) =∫ 50aexp(−s2 + a2)(4s)ds (2.40)= 2(1− exp(a2 − 502)). (2.41)20In this case, the reproductive value decreases with age and becomes 0 at maximum age.2.2.2 Multi Type ModelIn a population, individuals can be categorized into classes based on their physical features,habitat, and life stages. These classes are called types. The type of an individual changes withtime, for example in the case of immigration and emigration spatial location changes, in the caseof maturation life stage changes, in the case of birth events there are physical changes, and soon. These kinds of population models are multi-type or multi-state models. In this section weshow how to define reproductive value for multi-type models with time-independent transitionrates, using PF theory. Consider a population in which individuals are categorized in m types.Let n ∈ Z be the generation number. The densities of individuals in each type in generationn is denoted by vector X(n) = (x1(n), x2(n), x3(n), ..., xm(n)) in case of discrete-time modeldynamics. Suppose that Tij is the probability of transition from an individual of type i to typej in one generation given that it survives. Let T denote the m × m transition matrix. Forconsidering births, suppose bij is the number of offspring of type j produced by an individualof type i in one generation. Let B denote the m × m fertility matrix. Thus the density of anindividual in generation n+ 1 is,X(n+ 1) = X(n)(T +B) = X(n)A, (2.42)where coefficients of A are non-negative. The matrix A is known as the reproduction matrix.Type distribution in generation n can be written in terms of initial distribution as,X(n) = X(0)An. (2.43)Here, the transition rate and birth rate doesn’t dependent on the generation number, i.e., the ele-ments of matrixA remains constant with time. Hence, the reproductive value of the individualscan be determined by calculating the dominant non-negative eigenvalue λpf and correspondingnon-negative right eigenvector w, where wi denotes the relative reproductive value of an i typeindividual. Here, the non-negative left eigenvector v w.r.t λpf denotes the population distribu-21tion at equilibrium. We can estimate the unique left and right eigenvectors if and only if matrixA is a regular matrix. If matrix A is diagonalizable in terms of eigenvectors as basis vectors,we can write the solution of equation (2.42) as,X(n) =m∑i=0λni ciei, (2.44)where λi are the eigenvalues and ei are the corresponding eigenvectors. To analyze long-termdynamics, consider n → ∞ in equation (2.44). There are two cases for the eigenvalues, firstone is if the absolute value of all the eigenvalues is less than 1 then the population decays ex-ponentially. On the other hand, if the absolute value of some eigenvalues is greater than 1 thenthe population increases exponentially.Let’s consider a special case where the dominant eigenvalue, i.e., λpf = 1. In this case, thepopulation reaches an equilibrium state as n → ∞ and the equilibrium distribution is the lefteigenvector v of matrix A. This leads to another interpretation of the right eigenvector w (re-productive values), i.e., wi is the relative contribution of an individual of type i in the long-termpopulation.Example: An example of the multi-type model is a stage structured population model.Consider the life stages of an insect namely eggs, juveniles, and adults. The proportion of eggsthat get converted to juveniles is 0.6 in one generation and the proportion of juveniles whobecome adult is 0.2. An adult produces 5 eggs in one generation. The reproduction matrix forthis model is given by,A =0 0.6 00 0 0.25 0 0, where row and column are in the same order as the developmental stages. Let Xn be thepopulation size vector in nth generation. Suppose X0 = (10 50 20)T then the populationsize in nth generation is,Xn = X0An (2.45)22The dominant eigenvalue for matrix A is λpf = 0.84 and the reproductive value vector is,w =0.81.155 .In this case, only 20% juveniles can reach the adult stage which leads to less reproduction forinsects in the juvenile and egg stages. And an adult produces on average 5 eggs in its lifetime.Now, if we want to increase the insect population, we can perform sensitivity analysis on allthree transition rates to understand which life stage is of significance.2.3 Non-linear Population ModelsIn mathematical biology, a non-linear population model refers to a mathematical representa-tion of the dynamics of a population that takes into account interactions among individuals orinteractions with the environment. Predator-prey models [Berryman, 1992] are well-known ex-amples where non-linearity appears due to interaction among species. Broadly, non-linearityarises from density dependence, environmental conditions, and sexual reproduction.Density dependence is explored for plants in [Solbrig et al., 1988] and[Silva Matos et al., 1999] and for animals in [Longstaff, 1977],[Clutton-Brock et al., 1997],[Tanner, 1999] and [Cushing, 1988]. Competition for food, space, or other resources leads todensity dependence.Non-linearity due to environmental conditions can arise from many factors in-cluding fluctuations in natural resources [Henson and Cushing, 1997], climate change[Brodie and Post, 2010], seasonality, and weather. This type of non-linearity can be modeledvia density dependence because if the population is large it will impact the resource availabilitywhich in turn will impact the population density, which creates an environment feedback loop.The impact of considering sexual reproduction in the predator-prey model and host-parasitemodel is explored in [Doebeli and Koella, 1994].232.3.1 PeriodicityIn simple terms, when a set of events occur at fixed time intervals repeatedly then the event isperiodic in nature. In biological populations, periodicity is visible as breeding seasons, migra-tion patterns, and life cycles. During these periods vital rates and transition rates vary whichhas a significant impact on the growth of the population with time. In this section, populationmodels with periodicity in birth rate, death rate, and transition rate are considered for estimatingthe reproductive value [Jorgensen and Gromiec, 2016].Discrete Time ScenarioConsider a multi-type population in which individuals are categorized in m classes. The birthrate of an individual of type i in generation n is βi(n). Let B denote them×m fertility matrix.The probability for the transition of an individual from state i to state j in generation n is pij(n)given that it survives. Let T denote them×m transition matrix. Both the functions, birth rate,and transition rate, are periodic functions of n, i.e., within a period they are dependent onlyon the state of an individual. Suppose densities of individuals in generation n is denoted byX(n) = (x1(n), x2(n), ..., xm(n)). Assuming a discrete time setting, the density in generationn+ 1 is,X(n+ 1) = X(t)(T +B)(n) = X(n)A(n), (2.46)where coefficients of the matrix A are non-negative and function of the generation number andstate of an individual. Assuming that the period for elements of matrixA isN , i.e.,A(n+N) =A(n). Hence, starting from n = 0, the density in generation N is X(N) = X(0)(A(N −1)...A(2)A(1)A(0)) and if we go in multiples of period N then,X(nN) = X(0)(A(N − 1)...A(2)A(1)A(0))n (2.47)Assuming all the matrices (Ai) are irreducible and non-negative, then using the PF theorydefinition, the right eigenvector w of matrix (A(N − 1)...A(2)A(1)A(0)) corresponding to λpfis the reproductive value vector.24Continuous Time ScenarioConsider the setup of the age structured population model defined in section 2.1. Here, ratesare a function of time and age both, i.e., birth rate is β(a, t) and death rate is µ(a, t). The deathrate is defined as a periodic function of time with period d i.e µ(a + d, t + d) = µ(a, t) andbirth rate with the same period d i.e β(a + d, t + d) = β(a, t), where d ∈ Z. Suppose thepopulation density at time t of an age a individuals is u(a, t) then change in population densityis described by the equations,∂tu+ ∂au = −µ(a, t)u(a, t) (2.48)B(t) := u(0, t) =∫ A0β(a, t)u(a, t)da (2.49)u(a, 0) = φ(a) ∀a ∈ [0, A], (2.50)where B(t) is the total rate of births at time t. We can solve the above PDE using the methodof characteristics. Consider a cohort of population at time t0 and age a0 at time t0. Let u¯(h) =u(a0 + h, t0 + h) and µ¯(h) = µ(a0 + h, t0 + h). Then (2.48) is,du¯dh+ µ¯(h)u¯ = 0. (2.51)Solving the above ODE we get,u¯(h) = u¯(0) exp(−∫ h0µ¯(s)ds)(2.52)u(a, t) =φ(a− t) exp(− ∫ t0µ(a− t+ s, s)ds), a > tB(t− a) exp(− ∫ a0µ(s, t− a+ s)ds), a < t.(2.53)In this case, the distribution of trajectories from the renewal point is not identical as it dependson time, which is one of the criteria for defining the renewal state. Here, µ and β are dependenton the birth time of the individuals because individuals born at time t0 vs. individuals born attime s > t0 will have different life trajectories. Hence, we can’t use the renewal definition fordefining reproductive value.25In this particular model, birth age is a common feature for all individuals. Here, if try to definethe reproductive value in the same way as for time-homogeneous models as below then it mayfail to have the key property that it had in the case, namely, being invariant over time.h(a) =∫ ∞ap(s)p(a)β(s, t0 + a+ s)ds, (2.54)where p(s) = exp(− ∫ s0µ(x, t0 + x)dx) and t0 is birth time.2.3.2 Uniform Non-linearityIn some species, due to constraints on space and food resources, crowding or an increase inpopulation density impacts the growth of the entire population. Some of these factors impactall types of individuals in the same manner, also known as uniform dependence. In this section,uniform non-linearity in terms of dependence on total population size is discussed for continu-ous time and discrete time scenarios. One of the well-known examples is the population modelwith a logistic growth rate.Continuous time scenarioConsider an age structured population model with birth and death rate as a function of age andpopulation size at time t. Population size at time t is N(t) =∫ A0u(a, t)da, where u(a, t) ispopulation density for age a individuals at time t. Birth and death rates are defined as,β(a, t) = β(a) (2.55)µ(a, t) = µ(a) + f(N(t)), (2.56)where f : R 7→ R is a continuous function of total population size at time t. Note that if f is alinear function of N(t) then it is the logistic competition term [Guzowska et al., 2011] for thepopulation. The change in population density of this process is given by the following systemof equations,∂tu+ ∂au = −(µ(a) + f(N(t)))u(a, t) (2.57)26B(t) = u(0, t) =∫ A0β(a)u(a, t)da (2.58)u(a, 0) = φ(a) ∀a ∈ [0, A], (2.59)Note that, the right-hand side in equation (2.57) is the sum of µ(a)u(a, t) and f(N(t))u(a, t)i.e sum of the right-hand side of equation (2.9) for linear age structured model and the densitydependence term. Let u1(a, t) be the solution of equation (2.9), i.e., without non-linearity andu2(a, t) be the solution of equation (2.57), i.e., with non-linearity then,u2(a, t)u1(a, t)= exp(−∫ t0f(N(s))ds)(2.60)Here, also we can’t use the renewal theory definition because of time dependence. Hence, usingthe general definition reproductive value is defined as,h(a) =∫ Aap(x)p(a)β(x)dx, (2.61)where survival probability at age a is p(a) = exp(− ∫ a0µ(x) + f(N(t0 + x))dx) and t0 is thebirth time. The ratio of survival probability in the above equation can be written as,p(s)p(a)=exp(− ∫ s0µ(x) + f(N(t0 + x))dx)exp(− ∫ a0µ(x) + f(N(t0 + x))dx), s ≥ a (2.62)= exp(−∫ saµ(x) + f(N(t0 + s))dx). (2.63)Next, consider two particles v1 and v2 born at time t0. Suppose,N∗ is the equilibrium populationsize then the ratio of their reproductive value at age a for v1 and b for v2 such that b > a atequilibrium is defined as,hv1(a)hv2(b)=∫ Aaexp(− ∫ xaµ(s) + f(N∗)ds)β(x)dx∫ Abexp(− ∫ xbµ(s) + f(N∗)ds)β(x)dx(2.64)=∫ Aaexp(− ∫ xaµ(s)ds)β(x)dx∫ Abexp(− ∫ xbµ(s)ds)β(x)dx. (2.65)Note that, the ratio of reproductive value at equilibrium for a non-linear model is the same as27the one for linear age structured population models.Example: A typical example of uniform non-linearity is linear density dependence. Con-sider a population with birth rate as β(a) = 4a and death rate as µ(a) = 2a − N(t), whereN(t) is the population size at time t. Let the maximum age be 100 for all the individuals inthe population. Now, consider two particles born at time t = 0 namely v1 and v2. The ratio ofreproductive value of v1 at age 10 and v2 at age 20 at equilibrium is,hv1(10)hv2(20)=∫ 10010exp(− ∫ x102sds)4xdx∫ 10020exp(− ∫ x202sds)4xdx(2.66)=∫ 100104x exp(−x2 + 100)dx∫ 100204x exp(−x2 + 400)dx (2.67)= exp(−300)exp(102/2)− exp(1002/2)exp(202/2)− exp(1002/2) . (2.68)Discrete time scenarioConsider the multi-type population model from section 2.2 with the birth rate bij and transitionprobability pij as a function of population density at time t i.eN(t) = x1(t)+x2(t)+...+xm(t).Hence, the reproductionmatrix as a function of population density is denoted byA(N(t)). Now,population density vector at time t is X(t) = X(0)(A(N(t − 1))A(N(t − 2))...A(N(0))) =X(0)L(N), where X(0) is initial population vector. Assuming that there is a stable equilib-rium distribution v = limt→∞X(t) such that vL(v) = vλpf. The right eigenvector of matrixL(v) corresponding to λpf is the reproductive value vector for this model for the respectiveequilibrium state.In ecological population models, there is a possibility of multiple stable states, which re-sults in the variation of reproductive value. Early examples of multiple stable states are inmarine ecosystems, refer [Knowlton, 2004] for a detailed description of these models. Also, in[Diekmann et al., 2003], authors give a general framework to formulate and analyze the exis-tence of steady states for non-linear multi-type population models. Underlying this frameworkis the assumption that when environmental conditions remain constant, individuals are indepen-dent of one another. Moreover, interactions among individuals occur indirectly through theirimpact on the environment.282.3.3 State Dependent Non-linearityConsider the age structured model as defined in section 2.1 with birth and death rates as afunction of age a and population N(t) =∫ A0u(a, t)da at time t. Suppose population density attime t of age a individuals is u(a, t) then change in population density,∂tu+ ∂au = −µ(a,N(a, t))u(a, t) (2.69)u(0, t) =∫ A0β(a,N(a, t))u(a, t)da (2.70)u(a, 0) = φ(a) ∀a ∈ [0, A], (2.71)where N(a, t) = u(a, t)da i.e number of age a individuals. In this case, also, we can’t use therenewal theory definition because of time dependence. Hence, the expected number of offspringproduced by an individual in its remaining lifetime is given by the below formulation but it isnot an invariant reproductive value function as the birth rate of newborns is dependent on birthtime.h(a) =∫ Aap(x)p(a)β(x,N(x, t0 + x))dx, (2.72)where t0 is the birth time of an individual of age a. The probability to survive till age a forindividuals born at time t0 is defined as,p(a) = exp(−∫ a0µ(x,N(x, t0 + x))dx). (2.73)Here, the reproductive value and survival probability of an individual depends on age spe-cific population size. In some moth species, a high number of births will lead to less survivalrate for newborns because of limited resource availability, which is an example of state depen-dent non-linearity [Ruohomaki et al., 2003]. Numerical solutions of a non-linear age structuredmodel are discussed in [Cui and Chen, 2007] and the impact of density dependence in popula-tion dynamics of red deer is explored in detail by [Lande et al., 2006].29Chapter 3Stochastic Population ModelsThe previous chapter discusses deterministic population models, which enable us to deter-mine population size or distribution across states at any given time. In biological systems,the variability in demographic parameters and environmental factors leads to randomness inbirth events, death events, mutation, and migration patterns. By using a stochastic formula-tion, behavior at the individual level can be accommodated and analyzed for biological systemswith high variability in behavior across populations. Also, external factors such as environ-mental conditions, which have a significant and random impact on population models, can beincorporated into the stochastic formulation.Using PF theory in Chapter 2 for deterministic models, we get asymptotic population growthrate (λpf), reproductive value vector w, equilibrium population distribution v, and popula-tion size in large generation number n as a scaled form of equilibrium distribution v, i.e.,X(n) ∼ λnpfwTX(0)v. In stochastic models, the reproduction and survival processes followa distribution. For example, in deterministic models, we have a fixed birth rate whereas instochastic models we have an offspring probability distribution. Therefore, estimates of futurepopulation size given initial conditions vary, unlike deterministic models where we can predictexact future estimates. To describe long-term population dynamics in stochastic models, weanalyze and estimate the expected values and guess the most likely behavior (s) instead of anexact value.The reproduction matrixA from deterministic models becomes mean matrixM for stochas-30tic models, i.e., mij represents the expected number of offspring of type j produced by a typei individual. The key objects in stochastic models are the offspring distribution, the survivalprobability measure, and the distribution for the trajectory of an individual, i.e., the state ofan individual as a function of time. For example, let L be a random variable that representsthe number of offspring such that P (L = k) = pk for k ∈ N is the probability of having koffspring than the number of offspring of an individual in the deterministic case is equal toE[L]. For the stochastic models with mean matrixM and X(t) as population vector at time t,λpf is the asymptotic growth rate of E[X(t)], right eigenvector w is reproductive value and lefteigenvector v is asymptotic E[X(t)].In this chapter, we will prove the main result in several different contexts. The main resultis to show that for models that are critical in the sense that the expected population size neithergrows nor decays exponentially, the total reproductive value of the population is a martingale,which is a random variable that is constant on average in a sense that is made precise in Defini-tion 3.1.1. We will discuss the martingale property broadly in the context of total reproductivevalue and normalized population size. This chapter is organized as follows:• Branching processes - In section 3.1, we will discuss single type branching processesin continuous and discrete time scenarios and define their total reproductive value. Wewill also discuss the categorization of these processes w.r.t expected value of offspringdistribution. Further, theorems are given to understand the asymptotic distribution of thepopulation size.• Multi-type processes - In section 3.2, we will discuss the multi-type branching processes.We will define reproductive value using PF theory and discuss theorems to understandthe limiting behavior of the population.• Age structured populationmodel - In section 3.3, we will discuss a special case of the gen-eral branching process which is the age structured populationmodel. Wewill describe theconstruction of the process, the parameters required to analyze it, and its categorizationbased on the Malthusian growth rate. Further, we will show that total reproductive valueis a martingale only for critical processes and also prove the convergence in distribution31of a suitably rescaled population measure, conditioned on survival.Let’s start with understanding themathematical formulation of stochastic populationmodelsas branching processes.3.1 Branching ProcessesOne important type of stochastic process is the Markov process, in which the probability of anevent happening in the future depends only on the current state of the process instead of the entirehistory of the process. The branching process is a Markov process, in which each individualgives birth to a random number of offspring with respect to a probability distribution, i.e., theprobability of having i offspring at time t, where t denotes generation number in discrete timemodels and time in continuous time models. This probability distribution doesn’t depend on thehistory of the parent but it can vary with the state/type of the individual or time. In particular, ifoffspring probability distribution varies with time, it is called time in-homogeneous branchingprocess.3.1.1 Discrete Time ScenarioThe most common formulation of the branching process is the Galton-Watson (GW) process. Atypical discrete time GW process is defined with n as the generation number, L as the numberof offspring with P [L = k] = pk, i.e., probability of having k offspring. We assume thatindividuals are independent of each other and have identical offspring probability distribution.The expected number of offspring of an individual is denoted bym =∑∞k=0 kpk and populationsize at generation n as Zn. The expected population size in generation n is given by,Zn =Zn−1∑i=1Li(n) (3.1)E[Zn|Zn−1] =Zn−1∑i=1E[Li(n)] = mZn−1. (3.2)32Now, taking Z0 = 1 and by induction E[Zn|Z0 = 1] = mnZ0 = mn. The variance ofZn is given by the following formula (see Chapter 1 of [Athreya and Ney, 2004] for detailedcalculations) ,Var(Zn) =σ2mn−1(mn−1)m−1 ifm 6= 1,nσ2 ifm = 1,(3.3)where σ2 is the variance of offspring distribution. In the case of m = 1, variance in popula-tion size increases linearly with generation number, whereas when m < 1 then the variancedecreases exponentially with generation number, and when m > 1 then the variance increasesexponentially with generation number.Note that expected population size is dependent on m which implies long term dynamicsare dependent on the average number of offspring. Therefore, the GW process is classified onthe basis ofm in 3 categories:• sub-critical Process (m < 1): mean population size (E[Zn]) decreases with time,• critical Process (m = 1): mean population size (E[Zn]) remains constant with time,• super-critical Process (m > 1): mean population size (E[Zn]) increases with time.Note that the critical and sub-critical GW processes go extinct almost surely whereas the super-critical process may survive for a longer time [Athreya and Ney, 2004]. An individual is re-placed by its offspring (L) at the birth event and all the individuals are of the same type. Here,m is analogous to the PF root that is λpf, which gives the growth rate of the process. Since, themodel has only one type, the reproductive value vector w is one element vector with value 1.Before giving the next result, we give a precise definition of a martingale.Definition 3.1.1. [Durrett, 1999] Consider a probability space (Ω,F0, P ), where Ω is statespace, F0 is σ−field and P is the probability measure. Let Fn be a filtration, i.e., increas-ing sequence of σ−fields. A sequence of random variables Xn is said to be adapted to Fn ifXn ∈ Fn for all n. Here, If E[Xn+1|Fn] = Xn, then it is known as martingale.33Lemma 3.1.1. Let H(n) be the normalized population size in generation n, i.e., H(n) =m−nZn then H(n) is a martingale.Proof. Let initial population size be 1 which means H(0) = 1 then,E[H(n)|H(0) = 1] = E[m−nZn|H(0) = 1] = m−nE[Zn|Z0 = 1] (3.4)= m−nmn = 1 = H(0). (3.5)Hence, H(n) is a martingale. Note that, if m = 1, i.e., a critical branching process then popu-lation size Zn becomes a martingale.To better understand the long term population dynamics, we need to determine how fast thepopulation is decaying in sub-critical/critical cases and how fast the population is growing insupercritical case. These questions are answered by following theorems, detailed proof is givenin [Lyons et al., 1995].Theorem 3.1.1. [Kesten and Stigum, 1966] Suppose 1 < m < ∞ and H = limn→∞m−nZnthen, the following are equivalent:• P [H = 0] = q• E[H] = 1• E[L log+ L] <∞In supercritical GW process, ifE[L log+ L] <∞ then almost surelyH(n) has positive limitH that means Zn ∼ Hmn such that E[H] = 1. And if E[L log+ L] =∞ then limn→∞H(n) =0 almost surely. Hence, altogether, L log+ L is the decisive criterion here.Theorem 3.1.2. [Heathcote et al., 1967] The sequence {P [Zn > 0]/mn} is decreasing. Ifm < 1, then the following are equivalent:• limn→∞ P [Zn > 0]/mn > 0,• supE[Zn|Zn > 0] <∞,• E[L log+ L] <∞.34In the sub-critical GW process, survival probability decays exponentially with generationnumber n. Note that, even when the sub-critical process is conditioned on survival, it doesn’texplode as the expected value of Zn given Zn > 0 is bounded.Theorem 3.1.3. [Kesten et al., 1966]: Suppose that m = 1 and let σ2 := V ar(L) = E[L2]−1 ≤ ∞. Then we have:• Kolmogorov’s Estimate : limn→∞ nP [Zn > 0] = 2σ2• Yaglom’s limit law : If σ < ∞, then the conditional distribution of Znngiven Zn > 0converges as n→∞ to an exponential law with mean σ22. If σ =∞, then this conditionaldistribution converges to infinity.In the critical GW process, survival probability decreases as the process moves forward intime as denoted by Kolmogorov’s estimate. Note that, Yaglom’s limit law tells that in case of fi-nite variance for offspring distribution, population size, conditioned on survival, is proportionalto the generation number n.3.1.2 Continuous Time ScenarioConsider a branching process, in which birth events occur at random times, and the probabilityof having k offspring is denoted by pk with average m. The lifetime of each individual is anexponential random variable with the rate parameter β. We assume individuals are independent,have identical offspring probability distributions and an individual is replaced by its offspringat a birth event. As described in detail by [Athreya and Ney, 2004], let (Xt)t≥0 denote thecontinuous-time Markov chain, which represents this branching process. The transition matrixQ for the process is defined as,Q(0, 0) = 0, (3.6)Q(x, x+ k − 1) = βxpk, (3.7)Q(x, x) = −βx. (3.8)35Note that, Q is the transition rate matrix, i.e., Qij denote the jump rate from state i to state j. Inother words, Qij is the rate of change of transition probability from state i to state j.Next, we can define the expected value of change in population size in a small time intervalh > 0 as follows,E[Xt+h −Xt|Xt = x] = (0)(−βXth) +∞∑k=1(k − 1)pkβXth (3.9)=( ∞∑k=1kpk −∞∑k=1pk))βXth (3.10)= (m− 1)βXth. (3.11)Hence, we get an ODE for the expected value as,E[Xt+h −Xth∣∣∣Xt = x] = (m− 1)βXt (3.12)dE[Xt]dt= (m− 1)βE[Xt] (3.13)E[Xt] = E[X0] exp(β(m− 1)t). (3.14)Suppose X0 = 1 which implies that E[Xt|X0 = 1] = exp(β(m − 1))t. The variance ofXt is given by the following formula (see Chapter 5 of [Athreya and Ney, 2004] for detailedcalculations),Var(Xt) =σ2β−1(exp(2βt)− exp(βt)) if β 6= 0,σ2t if β = 0,(3.15)where σ2 is the variance of offspring distribution. Similar to the GW process, the variance inpopulation size is linearly dependent on the time when β = 0, and it is exponentially dependenton the time when β 6= 0.The continuous-time branching process is categorized into 3 categories as follows:• sub-critical process ifm < 1,• critical process ifm = 1 and36• super-critical process ifm > 1.Note that, the population decays exponentially in sub-critical cases whereas it grows exponen-tially in super-critical cases. In the critical case (m = 1), the expected population size is 1which doesn’t change with time.Lemma 3.1.2. Let H(t) = exp(−β(m − 1)t)Xt be the normalized population size at time t,then H(t) is a martingale.Proof. So, we want to show that E[H(t)|H(s)] = H(s) for some s < t. Using equation (3.14),E[H(t)|Fs] = E[exp(−β(m− 1)t)Xt|Fs] = exp(−β(m− 1)t)E[Xt|Fs] (3.16)= exp(−β(m− 1)t) exp(β(m− 1)(t− s))Xs (3.17)= exp(−β(m− 1)s)Xs = H(s) (3.18)In the critical case, population size Xt is a martingale.Example: (The general birth-death branching chain)Consider a continuous time branching process in which a particle gives birth to 2 or 0 childrenat the end of its lifetime. Hence, in each transition, the population either increases by 1 ordecreases by 1. Here, we are assuming that simultaneous birth events do not occur. Let α bethe death rate of each particle. Let p be the probability of producing 2 children and 1− p be theprobability of producing 0 children. Let Xt be the population size at time t and X0 = 1, thenthe transition rate matrix Q of the (discrete-time) embedded Markov chain is,Q(0, 0) = 1, Q(x, x+ 1) = p, Q(x, x− 1) = 1− p Q(x, x) = −αx. (3.19)Now, the expected number of offspring ism = 0∗(1−p)+2p = 2p and the expected populationsize is,E[Xt|X0 = 1] = exp(α(2p− 1)t). (3.20)37If p = 1/2, i.e., the probability of producing no children and 2 children is the same, then it is acritical process. Also, when p = 1, it is known as the Yule process [Durrett, 1999], and whenp = 0, it is a pure death chain [Durrett, 1999].3.2 Multi-type ProcessesIn this section, we discuss the branching processes for populations with multiple types of in-dividuals. Here, the type of an individual is its demographic and spatial features such as age,size, color, location, etc. In the biological population models, birth rate, death rate, and immi-gration parameters vary in accordance with demographic factors such as size, age, and genderof an individual. The dynamics of a population model depend on the type of individuals in thepopulation.Consider a population made up of individuals of m different types, in which each individ-ual is independent of each other, and offspring distribution of each type doesn’t change withgeneration. Let k = (k1, k2, ..., km)T be a vector-valued random variable that represents thenumber of offspring of an individual in each type, i.e., k1 is the number of offspring of type1. The offspring distribution for type i individuals is pi(k), where k is the offspring vector.Let mij be the expected number of offspring of type j produced by a type i individual in onegeneration, i.e., mij = E[Z1(j)|Z0(i) = 1] , where i, j ∈ {1, 2, ..., d}. Let M denote them×mmatrix withmij as elements. Suppose the population density vector in generation n− 1is Zn−1 = (Zn−1(1), ..., Zn−1(m))T , thenZn =m∑i=1Zn−1(i)∑j=1Lij(n− 1), (3.21)where Lij(n−1) is the vector valued random variable which represents the number of offspringof jth individual of type i in generation n − 1. Since all the individuals are independent andwithin each type, they follow the same offspring distribution, Lij(n− 1) follows the offspringdistribution of type i individuals. Let un = (un(1), un(2), ..., un(m))T denote the expectedpopulation size in generation n. Using the definition of Zn, expected population size of type `38individuals given Zn−1 is,E[Zn(`)|Zn−1] =m∑i=1Zn−1(i)∑j=1E[Lij,`(n− 1)] (3.22)=m∑i=1mi`Zn−1(i), (3.23)taking the expected value again we get the expected population size of type ` individuals as,E[Zn(`)] = E[E[Zn(`)|Zn−1]] =m∑i=1mi`E[Zn−1(i)] (3.24)un(`) =m∑i=1mi`un−1(i). (3.25)So, the expected population size vector in generation n is,un = un−1M (3.26)un = Z0Mn (3.27), where Z0 is the initial population vector. M is a non-negative matrix by definition ofmij andthe fact that Z > 0. Assuming M is also a regular matrix, the reproductive value is definedusing PF theory. The right eigenvector w is the reproductive value vector corresponding tothe dominant eigenvalue λpf. Multi-type branching process are classified into 3 categories inaccordance with λpf as follows:• if λpf < 1, then sub-critical,• if λpf = 1, then critical and• if λpf > 1, then super-critical.Note that the left eigenvector v w.r.t. λpf represents the asymptotic distribution of populationsize if the process survives.Lemma 3.2.1. The reproductive potentialH(n) = λ−npf Znw for a multi-type branching processis a martingale.39Proof. To show that H(n) is a martingale, we need to prove that E[Hn|Fn−1] = Hn−1, whereFn−1 is filtration for the process till generation n− 1. Using expected population size formulafrom equation (3.23) and eigenvector property of w,E[Hn|Hn−1] = E[λ−npf Znw|Fn−1] = λ−npf E[Zn|F ]w (3.28)= λ−npf Zn−1Mw = λ−npf Zn−1λpfw (3.29)= λ−(n−1)pf Zn−1w = Hn−1. (3.30)This implies that the expected value of reproductive potential doesn’t change with generationnumber, i.e., it is a martingale.Similar to the single type GW process, the following are the conditional limit theorems forthe multi-type branching process; see [Athreya and Ney, 2004] for detailed proof.Theorem 3.2.1. [Joffe and Spitzer, 1967] If λpf < 1 and E[kT log(k)] <∞ then,• limn→∞ P (Zn = y|Z0 = x, Zn > 0) = b(y), where b is a probability measure on statespace Rm and x is initial population vector and y is population vector in generation n.Here,m is the number of types in the population model.• limn→∞ λ−npf P (Zn > 0|Z0) = Q(0)ZT0 wTheorem 3.2.2. [Kurtz et al., 1997] There is a scalar random variableW such that for λpf > 1,limn→∞Znλnpf= Wv a.s.and P (W > 0) > 0 iffE[kT log(k)] <∞Theorem 3.2.3. [Joffe and Spitzer, 1967] If λpf = 1 andE[Z21 ] <∞, then ZTnw/n conditionedon Zn > 0, converges in distribution to the random variable with densityf(x) =1γ1exp(−x/γ1), x ≥ 040where,γ1 =vTwvTQ(w).Here,Q(w) is a multiplication of the left eigenvector v and diffusion coefficient of the offspringdistribution atw, i.e., the right eigenvector. Note that, in the case of single typeGWprocess, v =w = 1 with λpf = 1 and Q(w) = σ2/2, where σ2 is the variance of the offspring distribution.In words, after being conditioned on survival, reproductive potential in generation n scaleslike n multiplied with an exponential random variable. It is analogous to single type criticalGW process theorem 3.1.3.Example: Consider a basic Susceptible-Infected-Recovered (SIR) model used for studyinginfectious diseases. In this example, we are considering a discrete time setting. Let B be arandom variable that denotes the fraction of people infected by an infected person in unit timesuch that B follows a Poisson distribution with mean 0.02. Let H be a random variable thatrepresents the fraction of people recovered in unit time; it follows an exponential distributionwith themean 0.04. LetN be a random variable that denotes the fraction of people who becomessusceptible again; it follows an exponential distribution with a mean of 0.08. Therefore, themean matrix is,M =0.98 0.02 00 0.96 0.040.08 0 0.92 (3.31)Let Zn = (Zn(S), Zn(I), Zn(R)) be the population size vector in generation n. Suppose Z0 =(100, 4, 0) then,E[Zn] =[100 4 0]Mn.In this case, the dominant eigenvalue λpf = 1 and reproductive value vector w = (1, 1, 1) andasymptotic population distribution is v = (4, 2, 1).413.3 Age structured population modelIn this section, we will discuss the construction of the stochastic age structuredmodel and obtainthe associated martingale. Assuming that individuals are independent of each other with birthage 0, they follow identical reproduction and lifespan distribution. Hence, birth events are therenewal points of the process, and we can utilize the renewal theory definition to define repro-ductive value. We are using Ulam-Harris notation to formulate the age structured populationmodel as a branching process.Consider a set of finite sequences of positive integers as,Ω := {∅}⋃n∈N(N)n, (3.32)where N = {1, 2, ..} and ∅ denotes the root particles, i.e., the particle from which the processstarts. An element of Ω represents an individual in the process. For example, a particle u isrepresented as u = (∅, 3, 5, 6) implies that u is 6th child of 5th child of 3rd child of root particle.If a particle produces more than one offspring at a birth event, offspring are numbered from leftto right as shown in figure 3.1.Let state space for the process be S = [0,∞) and birth rate and death rate as a function ofage x ∈ S are β(x) and µ(x). For each particle u born at time s is given a triplet (ξu, Bu(t), σu),• a path ξu : [s,∞) → S, which gives information about the age of the particle at timet ≥ s,• a list Bu(t) which contains the time of birth events, label of offspring at birth and• σu, which is the lifetime of the particle.Let V (t) be the list of particles alive at time t. Let (X(t))t≥0 denotes the measure of the processat time t, i.e., X(t) =∑u∈V (t) δ(ξu(t)), where δ(.) is a δ measure, i.e., δ(x) is the measurewith,δ(x)(A) = 1(x ∈ A). (3.33)42Figure 3.1: Branching tree diagram.Let Nt be the population size at time t, i.e., Nt = |V (t)|. Let β(x) and µ(x) be age dependentbirth and death rates.The age structured branching process is categorized into 3 classes in accordance withMalthusian parameter α for the process as described in Section 2.2.1 of Chapter 2. Note that,expected growth α is 0 in critical cases, i.e., on average a particle replaces itself in its lifetime.It is analogous to m = 1 in a single type GW process and λpf = 1 in a multi-type branchingprocess.Here, reproductive value h(x) of a particle u of age x, which is also the eigenfunction oftransition operator in (2.14) and p(a) that is survival probability at age a are same as defined in(2.12) and (2.13) in Chapter 2.Lemma 3.3.1. Let H(t) be the reproductive potential at time t, which is defined asH(t) =∫ ∞0h(a)Xt(da),where Xt(da) is number of age da particles alive at time t then H(t) is a martingale if it is acritical process.43Proof. Note that, H(t) in lemma 2.2.2 is the same as E[H(t)] in the case of the stochastic agestructured population model. Now, according to lemma 2.2.2, E[H(t)] is constant over timefor the critical age structured population model, which implies thatE[H(t)|Fs] = E[H(s)|Fs] = H(s). (3.34)Hence, H(t) is a martingale if Xt is a critical process.Lemma 3.3.2. If (Xt)t≥0 is a critical age-structured branching process and Xt is countingmeasure, i.e., number of individuals alive at time t and letW be an exponential random variablewith mean 1. Then there exists a C > 0 such that,Xt([a,∞)])t→D CW∫ ∞ap(s)ds ∀a ∈ R+as t → ∞, where Xt([a,∞)) is the number of age ≥ a individuals at time t. Here p(s) isthe density function of the invariant measure of the population, which is analogous to the lefteigenvector in Perron-Frobenius theory.Proof. We will use theorem 2 from [Green, 1977], which states that conditioned on survivaldistribution of Xt/t converges in distribution to a exponential distribution as follows,Xtt→D exp( cχ¯), (3.35)where c > 0 is some constant that doesn’t depend on age andχ¯ =∫ ∞0E[χ(y)]dy. (3.36)In our case, let χx(t) be the characteristic function which 1 if particle x is alive and is of age≥ a else its value is 0.χx(t) =0, if σx < t,1, if ξx(t) ≥ a and σx > t,(3.37)44where σx is lifetime of the particle x. Using this definition of characteristic function expectedvalue of χ at age a is,E[χ(a)] = E[1(ξx(t) ≥ a and σx > t)] = p(a)1(ξx(t) ≥ a). (3.38)So,χ¯ =∫ ∞0E[χ(a)]da =∫ ∞ap(s)ds. (3.39)Using χ¯ in (3.35) with Xt conditioned on survival, we get the convergence to an exponentialdistribution as follows:Xt([a,∞))t→D exp( c∫∞ap(s)ds). (3.40)Now, rescaling to the exponential distribution with mean 1, it becomes,Xt([a,∞))t→D 1cexp(1)∫ ∞ap(s)ds. (3.41)45Chapter 4New characterizations of reproductiveValueIn Chapter 2, we have defined reproductive value in two equivalent ways: as the left eigenvectorcorresponding to the dominant eigenvalue using Perron-Frobenius theory and by counting thenumber of particles produced by an individual at renewal point using renewal theory. We alsodiscussed various deterministic population models and computed their reproductive values. InChapter 3, we discussed stochastic population models and their categorization in accordancewith the Malthusian growth rate. Further, we have shown that for critical population models,the total reproductive value is a martingale, which is computed by examining the mean behaviorof the stochastic model. We have also discussed the long-term dynamics of population modelsusing dominant eigenvalue and associated left and right eigenvectors, which represent equilib-rium distribution and reproductive value respectively. However, both of the above definitionsof reproductive value are limited to time-homogeneous population models, and the renewaldefinition requires a renewal point.In Chapter 3, we defined martingales for various branching processes. It was observed thatfor critical branching processes total reproductive value is a martingale. Specifically, for criticalstochastic population models, reproductive potential, i.e., total reproductive value at time t is46given by,Wt :=∫Sh(x)Xt(dx), (4.1)where S is the state space, h is the reproductive value function, andXt is the populationmeasureat time t, when Xt is a counting measure, then the above integral is equal to,∑v∈Vth(ξv(t)), (4.2)where Vt is the list of particles alive at time t and ξv(t) is the state of the particle v at time t.In the case of a time-inhomogeneous population model, we need a better characterizationor definition of reproductive value. Although in the case of population models with uniformnon-linearity, we defined the relative reproductive value at equilibrium which is the same as theone inhomogeneous case.In this chapter, we will obtain two characterizations of the ratio of reproductive values intime-homogeneous population models, which may be useful for defining the reproductive valuefunction for time-inhomogeneous population models. One is the relation with the limiting ratioof survival probability, and the other is the inverse ratio of martingale values with respect to thesize-biased measure, which is defined in section 4.2. In section 4.3, size-biased construction forthe continuous-time branching process and age-structured population models are discussed indetail. We will use reproductive potential as martingale to size bias the population model andto characterize the ratio of reproductive values of any pair of states x, y as the limiting inverseratio of the martingales corresponding to the process with a single particle initially at x and yrespectively, with respect to a specific coupling of those two processes (the notion of couplingis defined in section 4.2). The size-biased construction for the discrete-time branching processis discussed in detail by [Lyons et al., 1995].In this chapter, all the results are derived for critical stochastic processes. But, we knowthat a deterministic population model typically has a stochastic counterpart that we can use tocompute reproductive value, using these new characterizations.474.1 Reproductive value and survival probabilityConsider a single type critical GW process defined in section 3.1.1, i.e., the expected number ofoffspring is one (m = 1). In this case, the reproductive value in generation n is the populationsize in generation n.Lemma 4.1.1. Let (Zn)n≥0 be a critical single-type GW process. Consider another copy of theprocess (Sn)n≥0 such that Z0 = q1 and S0 = q2, where q1.q2 ∈ N \ {0}. Suppose H(n) isreproductive potential in generation n for Sn process and V (n) be the reproductive potentialfor Zn process then,V (0)H(0)=q1q2, (4.3)andlimn→∞nP (Zn > 0|Z0 = q1)nP (Sn > 0|S0 = q2) =V (0)H(0)(4.4)Proof. Let sn(z) be the survival probability given that initial population size is z, i.e., sn(z) =P (Zn > 0|Z0 = z) then extinction probability is ρn(z) = 1 − sn(z). Since all individualsare independent, so if the process starts with > 1 individuals, then ρn(z) = ρn(1)z. Therefore,sn(z) can be written as,sn(z) = 1− ρn(z) = 1− ρn(1)z = 1− (1− sn(1))z. (4.5)Now, from theorem 3.1.3, we know that,limn→∞nsn(1) =2σ2, (4.6)where σ2 is the variance of offspring distribution. Since the survival probability of a criticalGW process is inversely proportional to the generation number, which means sn(1) < 1 for48large n, and more importantly sn(1)→ 0 as n→∞. Therefore, (4.5) is approximated as,sn(z) ∼ 1− (1− zsn(1)) = zsn(1) (4.7)Hence, the ratio of survival probabilities with different initial numbers of individuals is,limn→∞nP (Zn > 0|Z0 = q1)nP (Sn > 0|S0 = q2) =q1q2, (4.8)which is the ratio of reproductive values.Next is the generalized version of lemma 4.1.1 for the age-structured branching process.Theorem 4.1.1. Let (X(t))t≥0 be an age-structured critical branching process such that h(x)is the reproductive value of age x particle. Suppose there are two lineages; one starts with agex particle and the other start with age y particle. Let B(t) be the number of birth events tilltime t and let τi denote the birth time for ith offspring of the root particle of a process X(t). If∃α > 0 such that µ(a) ≥ α for all a ∈ R+, E[Bx(t)] <∞ , E[By(t)] <∞ and τi is boundedthen,limt→∞Sx(t)Sy(t)=h(x)h(y), (4.9)where Sx(t) and Sy(t) is survival probability of lineage starting with age x particle and age yparticle respectively.Proof. Let ρx(t) denote the probability of extinction by time t for a population starting with anindividual v of age x then Sx(t) is given by,Sx(t) = 1− ρx(t). (4.10)The lineage of v is extinct by time t iff the particle v, and the lineages of offspring of v born bytime t, are dead by time t. Let B(t) denote the number of offspring of v in [0, t] then,ρx(t) = P (ζv ≤ t and (ζ{vi} ≤ t− τi) ∀i ≤ B(t)), (4.11)49where ζv is the death age for particle v, ζ{vi} is the lifespan of the lineage of ith offspring and τiis birth time of ith offspring. Here, τi’s are labeled in increasing order. It is helpful to re-writeρx(t) using indicator functions as follows:ρx(t) = E[1{ζv ≤ t}ΠB(t)i=1 1{ζ{vi} ≤ t− τi}] (4.12)Let Ft be the σ algebra up to time t generated by the following random variables: lifespanof particle v till time t and offspring measure of particle v up to time t, i.e., birth time of ithoffspring of particle v, i.e., τi ≤ t. Here, particles evolve independently, and the birth eventis the renewal point implying that the life trajectory of the offspring after birth is independentof what happened before the birth event. Therefore, ζv and τi ∈ [0, t] are Ft measurable whileζ{vi} is independent of Ft for all i ≤ B(t). Now, we can write ρx(t) as,ρx(t) = E[E[1{ζv ≤ t}ΠB(t)i=1 1{ζ{vi} ≤ t− τi}|Ft]] (4.13)= E[1{ζv ≤ t}E[ΠB(t)i=1 1{ζ{vi} ≤ t− τi}|Ft]] (4.14)= E[(1− 1{ζv > t})ΠB(t)i=1 P (ζ{vi} ≤ t− τi |Ft)] (4.15)= E[ΠB(t)i=1 P (ζ{vi} ≤ t− τi |Ft)]− E[1{ζv > t}ΠB(t)i=1 P (ζ{vi} ≤ t− τi |Ft)]. (4.16)Since, ζ{vi} is independent of Ft we can write the extinction probability of the lineage of ithoffspring of v as ρ0(t − τi) = P (ζ{vi} ≤ t − τi |Ft). Let ρ∗x(t) = E[ΠB(t)i=1 ρ0(t − τi)] andRx(t) = E[1{ζv > t}ΠB(t)i=1 ρ0(t− τi)] then expression for ρx(t) is,ρx(t) = ρ∗x(t)−Rx(t) (4.17)Note that, 0 ≤ Rx(t) ≤ P (ξv > t) which follows from the observation that ∀i, 0 ≤ ρ0(t −τi) ≤ 1, which means 0 ≤ ΠB(t)i=1 ρ0(t − τi) ≤ 1. Hence, Rx(t) is at most P (ξv > t). Notethat survival probability for a single particle, i.e., v decreases exponentially as P (ξv > t) =exp(− ∫ t0µ(v(s))ds), where µ(s) is the death rate at age s. Assuming µ(a) > α for all a ∈ R+implies that 0 ≤ Rx(t) ≤ P (ξv > t) ≤ exp(−αt). For a general critical branching processstarting with a particle of age x, ρx(t) ∼ c/t as t → ∞ for some c that depends on offspring50distribution as shown in [Green, 1977]. Hence, Rx(t) ≤ exp(−αt) = o(Sx(t)).Next, take ρ∗x(t) and using the fact that S0(t− τi) = 1− ρ0(t− τi), we getρ∗x(t) = E[ΠB(t)i=1 (1− S0(t− τi))] (4.18)Lets consider the expansion of product inside expectation with ai = S0(t − τi) for i ∈{1, 2, ..., B(t)}, which is a inclusion exclusion series,ΠB(t)i=1 (1− ai) = 1−B(t)∑i=1ai +∑1≤i<j≤B(t)aiaj −∑1≤i<j<k≤B(t)aiajak + ... (4.19)Now, looking at the third term in the expansion,∑1≤i<j≤B(t)aiaj =12B(t)∑i=1aiB(t)∑j=1aj − 12B(t)∑i=1a2i (4.20)=12( B(t)∑i=1ai)2− 12B(t)∑i=1a2i (4.21)Since, 0 ≤ ai ≤ 1 for all i ∈ {1, 2, .., B(t)} and (4.19) is a alternating series,0 ≤ ΠB(t)i=1 (1− ai)−(1−B(t)∑i=1ai)≤ 12( B(t)∑i=1ai)2− 12B(t)∑i=1a2i (4.22)<12( B(t)∑i=1ai)2. (4.23)Next, using Jensen’s inequality we can say that,( 1B(t)B(t)∑i=1ai)2≤ 1B(t)B(t)∑i=1a2i (4.24)( B(t)∑i=1ai)2≤ B(t)B(t)∑i=1a2i < B(t)B(t)∑i=1ai. (4.25)51Using ai = S0(t− τi),B(t)∑i=1ai =B(t)∑i=1S0(t− τi) (4.26)=B(t)∑i=1S0(t− τi)t− τi (t− τi). (4.27)Using the fact that tS(t) ∼ c as t→∞ for some c as shown in [Green, 1977] for critical generalbranching process, i.e., S0(t− τi)(t− τi)→ c as t→∞. Assuming that τi is bounded, i.e., anindividual doesn’t reproduce late in its lifetime, which means that limt→∞ 1/(t − τi) → 0 andalso assuming that E[B(t)] <∞ we get,limt→∞B(t)∑i=1ai → 0. (4.28)Hence,(∑B(t)i=1 ai)2= o(∑B(t)i=1 ai) which implies that as t→∞,ΠB(t)i=1 (1− ai) ∼ 1−B(t)∑i=1ai. (4.29)Hence, approximation of ρ∗x(t) is,ρ∗x(t) ∼ E[1−B(t)∑i=1S0(t− τi)](4.30)= 1− E[ B(t)∑i=1S0(t− τi)](4.31)= 1− E[S0(t)B(t)∑i=1S0(t− τi)S0(t)]. (4.32)Now, collecting the terms for Sx(t) assuming S0(t) is constant,Sx(t) = 1− (ρ∗x(t)−Rx(t)) (4.33)Sx(t) ∼ S0(t)E[ B(t)∑i=1S0(t− τi)S0(t)]+Rx(t). (4.34)52Note that for large t, the ratio inside expectation that is S0(t−τi)/S0(t) denotes the survivalprobability of a lineage started at time τi given that it survives till time t. Let’s look at the ratioof survival probabilities,S0(t− τi)S0(t)=S0(t− τi)S0(t)t− τittt− τi (4.35)=(t− τi)S0(t− τi)tS0(t)tt− τi . (4.36)If the variance of offspring distribution is finite for a critical general branching process, thentS(t) ∼ c as t → ∞ for some c which depends on offspring distribution as shown in[Green, 1977]. So, using tS(t) ∼ c and noting that t/(t − τi) ∼ 1 as t → ∞ in (4.36) weget the pointwise convergence as below,S0(t− τi)S0(t)∼ 1. (4.37)Again, assuming that τi’s are bounded, we can strengthen this result to convergence in expec-tation. Therefore first term in (4.34) ,limt→∞E[ B(t)∑i=1S0(t− τi)S0(t)]= Ex[B(∞)] (4.38)Hence,limt→∞Sx(t)S0(t)= Ex[B(∞)]. (4.39)Similarly for other lineage starting with particle of age y, limt→∞Sy(t)S0(t)= Ey[B(∞)]. Sincethe expected number of birth events for an individual in its remaining lifetime is equal to thereproductive value.limt→∞Sx(t)Sy(t)=Ex[B(∞)]Ey[B(∞)] =h(x)h(y). (4.40)53We know that the Malthusian growth parameter α is greater than 0 and less than 0 for super-critical and sub-critical processes (Chapter 3). This parameter is calculated using the belowrelation,∫ ∞0exp(−αt)p(t)β(t)dt = 1, (4.41)where p(t)β(t) is the net birth rate at age t. Let I(0) be the expected number of offspring of aparticle in its entire lifetime which is given by,I(0) =∫ ∞0p(t)β(t)dt. (4.42)In the case of super-critical branching process I(0) > 1, which implies that α > 0. Now, if weincrease the number of death events by modifying p(t) uniformly across all individuals in thepopulation such that I(0) = 1, then (4.40) will be valid for the super-critical process also.Similarly, in the case of a sub-critical branching process I(0) < 1, which implies thatα < 0. Here, if we increase the number of birth events by changing the β(t) uniformly acrossthe population such that I(0) = 1, then (4.40) will be valid for the sub-critical process also.4.2 Size BiasingConsider a population with N number of families with different family sizes (number of chil-dren). For estimating the average of number of children per household, it is decided to randomlyselect n children. Then data for number of siblings of the selected children from sample is col-lected to estimate the average family size. This average family size represents inflated versionof actual family size, i.e., it is size biased. Since, the probability of selecting a child from alarge family is higher as compared to small families, it will result in overestimation.Definition 4.2.1. [Arratia et al., 2019] Let X be a non-negative random variable with E[X] ∈(0,∞) and X∗ be the random variable with size biased distribution of X . If X is a discrete54random variable, thenP (X∗ = x)P (X = x)=xE[X], (4.43)and if X is a continuous random variable, thenP (X∗ ∈ dx)P (X ∈ dx) =xE[X]. (4.44)The characterization for probability mass function f(x) in case of discrete random variable Xor probability density function f(x) in case of continuous random variable X is given by,fX∗(x) = f(x)xE[X]. (4.45)Let’s see how size biasing appears while estimating the average family size. Suppose thereare n different family sizes in the population such that the total population size is N . Let Xidenote the size of ith family for i ∈ {1, 2, ..., n} such that Xi’s are i.i.d random variables andS =∑ni=1Xi. Let I denote the index for the family size of the selected sample such thatY := XI then,P (I = i) =XiN. (4.46)Now, we want to know the probability that selected sample Y belongs to the family of size I ,i.e.,P (YI = x) =n∑i=1P (Xi = x, I = i) (4.47)=n∑i=1P (I = i|Xi = x)P (Xi = x) (4.48)=n∑i=1E[ xx+ S ′]P (Xi = x), (4.49)where S ′ = S −Xi = S − x and S ′ is sum of n− 1 i.i.d random variables. Now, we can writeS ′/(n − 1) → E[X] as n → ∞. If Xi = x is fixed then n/(x + S ′) ∼ n/E[X] as n → ∞ in55probability. Here, our sample space is discrete, and convergence in expected value holds for afixed x, then it holds after averaging over x. Hence, for large n,P (YI = x) =n∑i=1xnE[X]P (Xi = x) (4.50)= nxnE[X]P (Xi = x) =xE[X]P (Xi = x). (4.51)In population models, size-biased birth rates are used to construct the size-biased version ofa population process. This construction provides a convenient framework for characterizingthe long term behavior of the process in terms of survival probability, which we will discussin section (4.2.2). Since there is a relation between ordinary and size-biased distribution, byanalyzing the population process in the size-biased version, conclusions can be drawn for theoriginal process.4.2.1 Martingale change of measureDefinition 4.2.1 of size biasing concerns a single random variable. In this section, we willunderstand how a population process is size biased, which consists of a collection of randomvariables.Consider a critical branching process (Xt)t≥0 on measurable space (Ω,F , P ), where Ω isZ+, F is σ-algebra of the process and P is probability measure. Here,Xt is the population sizeat time t. Let Wt denote the total reproductive value of the critical process (Xt)t≥0 at time t,which is a non-negative martingale. Using this martingale, we can define a new measureQt onthe measurable space (Ω,Ft) as follows,Qt(A) =∫AWt(w)W0(w)dPt(w) for A ∈ Ft, (4.52)whereQt and Pt represents measureQ and P restricted to Ft and Ft is σ-algebra till time t. Atfixed time t, measure Qt is absolutely continuous w.r.t Pt, i.e., Q(.|Ft) << P (.|Ft) becauseits Radon-Nikodym derivative is Wt/W0 which has E[|Wt/W0|] = E[Wt/W0] < ∞ by non-negativity ofW and the martingale property [Durrett, 2019].56The martingale property ofWt implies that (Wt)t≥0 is adapted to a filtration (Ft)t≥0, i.e.,EP [Wt|Fs] = Ws for t > s, (4.53)where EP is expected value w.r.t P measure such thatWt → W∞ P almost surely as t→∞.Since, Wt is a martingale with respect to Ft under measure P , the measures Qt are consistentwith respect to Ft as shown below,Qt(A|Fs) =∫AWt(w)W0(w)dPt(w|Fs) (4.54)=∫AWs(w)W0(w)dPs(w) (4.55)= Qs(A), (4.56)where second equality is due to the property of conditional expectation as stated below.Theorem4.2.1. [Durrett, 2019] Consider a probability space (Ω,F0, P ), a sigma fieldF ⊂ F0and a random variable X measurable with respect to F0 with E[X] < ∞. If F ⊂ G andE[X|G] ∈ F then E[X|F ] = E[X|G].In our case, Fs ⊂ Ft as t > s and random variable isWt(w)/W0(w), which is a martingale(4.53) w.r.t measure P .Next, by Carathéodory’s extension theorem, there exists a unique probability measureQ onF such that Q|Ft = Qt. Here notation P |Ft and Q|Ft represents restriction of measure P andQ to σ-algebra Ft. The probability measure Q is known as the martingale change of P .If we have a non-negative martingale for the underlying branching process, we can define anew probability measure using equation (4.52). From Chapter 3, we know that the total repro-ductive value of various branching processes is a martingale. Now, we can use that to definea new measure. Before that, let’s define the probability measure for the ordinary branchingprocess. Consider a single type GW process from (Zn)n≥0 Chapter 3 with Z0 = 1. Let Ln,i bea random variable that denotes the number of offspring of ith individuals in generation n, andit follows the offspring distribution pk. All individuals in the population are independent and57follow the offspring distribution pk. If Zn is the population size in generation n then,Zn+1 =Zn∑i=1Ln,i. (4.57)Now, for this process probability measure P is,P (Zn+1 = j|Zn = i) =p∗ij if i ≥ 1, j ≥ 00 if i = 0, j ≥ 0,(4.58)where p∗ij is i-fold convolution of pk for k ∈ N.Note that, H(n) = m−nZn is a martingale under P measure as shown in lemma 3.1.1. Inthis case, if P is the probability measure for the GW process then the probability measure forsize biased GW process will be,dQdP∣∣∣Fn=ZnE[Zn]=Znmn. (4.59)In a similar way, for multi-type branching process (Zn)n≥0 from chapter 3 with M as meanmatrix and w as reproductive value vector, size biased measure is defined as,dQdP∣∣∣Fn=ZTnwZT0 w, (4.60)where the right-hand side is a martingale under measure P as shown in lemma 3.2.1. Similarly,for the continuous time branching process (Xt)t≥0 with β as birth rate and m as the expectednumber of offspring, the probability measure for size biased version is given by,dQdP∣∣∣Ft=XtE[Xt]=Xtexp(β(m− 1)t) , (4.61)where the right-hand side is a martingale under measure P (lemma 3.1.2).584.2.2 Spine process constructionIn earlier studies [Feller, 1991] and [Athreya and Ney, 2004], branching processes are studiedusing analytical methods, which are not very transparent. In recent years, intuitive probabilisticmethods are getting used to study the dynamics of the branching process. One such method,known as the spinal decomposition, has seen extensive use in recent work on branching particlemodels [Boenkost et al., 2022],[Foutel-Rodier and Schertzer, 2022], [Lyons et al., 1995]. Thespinal decomposition is a probability measure on the set of spacetime trees with a distinguishedpath, and on the set of (unlabelled) trees, it projects to the size-biased measure. In this method,the branching process is defined in a manner such that there is a trajectory that never endsthat lives forever. This trajectory is known as the spine of the process. In the case of criticaland sub-critical branching processes, the spine is the unique infinite lineage (spacetime path ofdescendants). Using a martingale change of measure to define a new measure, the transformedprocess can be decomposed in probability measure as an immigration process along this spinepath. This reduces the study of the branching process from many random trajectories to onerandom trajectory. The behavior of spine trajectory is utilized to study the long-term dynamicsof the branching process.Consider a continuous-time branching process (Xt)t≥0 in which the lifetime of an indi-vidual is an exponential random variable with rate β and pk as a probability distribution foroffspring, i.e., P (L = k) = pk, where L is the random variable denoting number of offspringsuch that E[L] = m. All individuals are independent, and an individual is replaced by theiroffspring at the time of death. Since more branches are generated at a death event, β is alsoknown as the branching rate. Let (X∗t )t≥0 denote the size-biased version of this branching pro-cess. The steps to construct a size-biased branching process denoted by (X∗t )t≥0 are as follows[Hardy and Harris, 2006] :1. the process starts with an ancestor or a root particle labeled as∅with size-biased branch-ing ratemβ2. root particles produced k offspring with size-biased offspring distribution. Let L∗ be the59size-biased version of L, thenP (L∗ = k) =kpkm(4.62)3. one offspring out of k is selected uniformly at random, which will follow size-biasedoffspring distribution and size-biased branching rate, and offspring distribution for re-maining k − 1 offspring is the same as of L, i.e., pk. The chosen individual is marked asv04. repeat step 2,3 and mark the chosen individuals as v1,v2 and so on.We are biasing the offspring distribution and branching rate only for the chosen individual. Thelist of chosen individuals (v0, v1, ...) forms a path that never ends because P (L∗ = 0) = 0. Thisdistinguished path is known as the spine of the process.Let P denotes the probability measure of (Xt)t≥0 and Q∗ be the probability measure of(X∗t )t≥0. Let random variable nt denote the number of birth events on the trajectory of a particlein ordinary process (Xt)t≥0, and it follows Poisson distribution with rate β. Let random variablen∗t be the number of birth events on the spine trajectory in size-biased process (X∗t )t≥0, and itfollows size-biased Poisson distribution with rate mβ. Suppose there are n birth events on thespine path till time t then,P (n∗t = n) = exp(−mβt)(mβt)n1n!,and (4.63)P (nt = n) = exp(−βt)(βt)n 1n!. (4.64)Hence, the extra factor for branching on the spine path is exp(−(m− 1)βt)mn.Next, measure Q∗ for the size-biased process with a distinguished spine path is defined as,Q∗(., (v0, v1, ..., vnt−1)) = exp(−(m− 1)βt)mntΠv∈(v0,...,vnt−1)( km1k)P(.). (4.65)Note that, the difference between the two processes is because of accelerated branching rate andsize-biased offspring distribution, as other individuals grow in the same way as the unbiased60process. Therefore, the extra terms required are:• exp(−(m− 1)βt)mnt to accommodate the accelerated branching rate and• kpk/m ∗ 1/k = 1/m to account for bias in offspring distribution, and selecting the spineparticle at each birth event on the spine path. For nt birth events, it will be 1/mntTherefore,dQ∗dP∣∣∣Ft= exp(−(m− 1)βt)1((v0, ..., vnt−1)) (4.66)Note that, the right-hand side in the above relation is the inverse of the expected population sizefor a continuous time branching process with a labeled spine path (section 3.1.2).Next, consider the case when we know that process is size biased, but we don’t have infor-mation about the spine path, i.e., (v0, v1, ..., vnt−1). In other words, the spine path is not labeled.So, at time t, all Xt individuals are equally likely to be the spine particle. Let Q be the prob-ability measure for this process in which the spine path exists but is not labeled. It is definedas,Q(.) =∑v∈XtQ∗(., (v0, v1, ..., vnt−1)) (4.67)= exp(−(m− 1)βt)XtP(.). (4.68)Hence, we get the relation (4.61) between Q and P probability measures as,dQdP∣∣∣Ft= exp(−(m− 1)βt)Xt. (4.69)4.2.3 Age Structured Population ModelIn this section, we will define the spinal representation for the age structured branching process.To define the spinal representation, we need to construct a size-biased process using steps fromthe previous section and label the spine particles at each birth event on the spine trajectory ofthe process.61Consider the age structured process defined in section 3.3 and let L be the random variabledenoting number of offspring. We only consider the case where an individual produces only1 offspring with probability 1, i.e., P (L = 1) = 1. To construct the spinal representation, weneed to size bias two parameters as follows,• size biased branching rate, which ism = 1 and• size-biased offspring distribution, which is 2q1/(m+ 1) = q1, where q1 = 1.The construction of a size-biased age structured population model starts with a root particle ofage x ∈ R+ at time s. At the first birth event of root particle at time τ1, it produces 1 offspring.Then the label of the spine particle is randomly passed on to the newborn or remains with theparent, with a probability proportional to their reproductive values. Let Bv(t) be the number ofbirth events till time t on the spine path. Let v = (v0, v1, ..., vB(t)−1) be the list of spine particlestill time t.Let P be the probability measure for the unbiased process, then probability measure Q∗ forthe size-biased process with distinguished spine path is given by,Q∗(., v) = (2)B(t)Πi∈{0,1,...,B(t)−1}(22h(vi)h(0) + h(vi−1))P (.) (4.70)= Πv∈v( 2h(ξv)h(0) + h(f))P (.), (4.71)where f is the age of the parent spine particle at respective birth event.Taking the marginal distribution of Q∗ w.r.t all possible v. Suppose set of list of spineparticles till time t is denoted by G(t) then,Q(.) =∑v∈G(t)Πv∈v( 2h(ξv)h(0) + h(f))P (.) (4.72)=1h(v0)∑u∈Xth(ξu(t))P (.). (4.73)Note that, 1h(v0)∑u∈Xt h(ξu(t)) is a martingale. It is denoted byWx(s, t), which means martin-62gale of a process till time t, which starts with a particle of age x at time s. Hence,dQxdPx∣∣∣Ft=Wx(t)Wx(0), (4.74)whereWx(t) := Wx(0, t).Here, spine movement is deterministic for the part when the spine particle is aging withtime, and then it jumps to age 0 at time t with rate λ(t) given by,λ(t) =h(0)h(ξv(t)) + h(0)h(ξv(t)) + h(0)h(ξv(t))β(t) (4.75)=h(0)h(ξv(t))β(t), (4.76)where ξv(t) is the age of parent spine particle and β(t) is ordinary birth rate at time t.The stationary distribution pi for the spine process starting with a particle of age 0 ispi(t) = p(t)pi(0) = exp(−∫ t0λ(s)ds)pi(0). (4.77)Lemma 4.2.1. For an age structured population model, the stationary distribution pi of a spineprocess starting with a particle of age 0 at t = 0 is proportional to the product of eigenfunctionh and eigen probability measure p which are defined as follows,h(a) =∫ ∞ap(s)p(a)β(s)ds (4.78)p(a) = exp(−∫ a0µ(s)ds) (4.79)Proof. We know that pi(t) = exp(− ∫ t0λ(s)ds). Let’s differentiate log(pi) w.r.t t.log(pi(t))′ =1pi(t)pi(t)′ =−1pi(t)λ(t)pi(t) (4.80)= −λ(t) = − h(0)h(ξv(t))β(t) (4.81)= −β(t)h(t)(4.82)63The last equality is because the population model is a critical branching process which meansh(0) = 1. Also, assuming that the spine process starts with a particle of age 0 at t = 0, its ageat time t is t, i.e., ξv(t) = t. Next, the product of eigenfunction and eigen probability measureis,h(t)p(t) =∫ ∞tp(s)β(s)ds. (4.83)Differentiating log(h(t)p(t)) w.r.t t,log(h(t)p(t))′ = − 1h(t)p(t)p(t)β(t) =−β(t)h(t). (4.84)Hence log(pi(t))′ = log(h(t)p(t))′. So, if initial conditions at t = 0 are the same then thestationary distribution of the spine is equal to h(t)p(t).From equation (4.82) and (4.84),pi(t) = c1 exp(−∫ t0β(s)/h(s)ds) (4.85)h(t)p(t) = c2 exp(−∫ t0β(s)/h(s)ds), (4.86)where c1, c2 > 0 are some constants. Hence, if c1 = c2 = c > 0 then pi(t) is equal to h(t)p(t)up to a positive constant c.Now, lets understand how does reproductive value behaves w.r.t β and p, starting withdifferentiating both sides of the equation (4.83),h(t)p′(t) + h′(t)p(t) = −p(t)β(t) (4.87)h′(t) =1p(t)(−p(t)β(t)− h(t)p′(t)) = −β(t) + µ(t)h(t) (4.88)h′(t) = µ(t)(h(t)− β(t)µ(t)). (4.89)Here, we have an ODE for h(t), and we know that for critical process h(0) = 1, which is aboundary condition for the above ODE, and the initial condition is h(∞) = 0. Let’s consider a64special case when β(t) = µ(t) = 1 for all t ∈ R+, which is an equilibrium point. So, the ODEsystem becomes,h′(t) = h(t)− 1,h(0) = 1h(∞) = 0.(4.90)The solution is h(t) = 1 for all t ∈ R+.4.3 CouplingCoupling is a technique in probability theory to compare two random variables X and Y bycreating a random vectorW such that its marginal distributions corresponds to the distributionsof X and Y .Definition 4.3.1. A coupling of two random variables X and Y on the same measurable space(Ω,Σ) is any pair of random variables (Xˆ, Yˆ ) on the product measure space (Ω × Ω,Σ⊗Σ)whose marginals have the same distribution as X and Y ,Xˆ =D X, Yˆ =D Y.Note that, the probability distribution Pˆ of (Xˆ, Yˆ ) is a coupling of probability distributions ofX and Y that are P and P ′ respectively.Example 1 (A game of two biased coins): Suppose you have two coins, and the outcomeof tossing them is X and Y respectively. The outcome head is denoted by 1, and the tail isdenoted by 0. The probability of getting heads for one coin P (X = 1) is p, and for another coinP (Y = 1) is q such that q > p. The probability of getting heads on the second coin is higher thanon the first coin. But the question arises, is there a way to throw two coins such that outcome isalways ordered, i.e., X ≤ Y ? Yes, it is possible by coupling the probability distributions of X65and Y . To couple these random variables, a third coin with an outcome denoted by Z is createdsuch that P (Z = 1) = q−p1−p = q′. Now, toss a first and third coin and assign the maximum oftheir outcomes as Y ′ = max(X,Z). Note that the probability distribution of Y ′ is the same asthe distribution of Y .P (Y ′ = 1) = P (X = 1, Z = 0) + P (X = 0, Z = 1) (4.91)= p(1− q′) + (1− p)q′ (4.92)= q (4.93)Hence, the coupling (X,Y ′) ensures that X ≤ Y ′ for all tosses with marginal distributions ofX and Y respectively.Another illustration is an application of coupling for proving that the asymptotic probabilitydistribution of two random walks is the same irrespective of their starting position in the 1dimension, where the difference in starting positions is an even number. Before that, let’sdefine total variation distance and understand its relation with coupling.Definition 4.3.2. The total variation distance (TV) between two probability distributions P andP ′ on measurable space (Ω,Σ) is,||P (X = x)− P ′(Y = x)||tv = maxx⊆Ω|P (X = x)− P ′(Y = x)|Total variation distance represents how close the distribution of two random variables Xand Y on the same probability space. The relation of TV with coupling is,||P (X = x)− P ′(Y = x)||tv ≤ inf{Pˆ (Xˆ 6= Yˆ )}, where (Xˆ, Yˆ ) is coupling of the random variables X and Y . Another way to characterize TVis,||P (X = x)− P ′(Y = x)||tv = 12∑x∈Ω|P (X = x)− P ′(Y = x)|An important coupling inequality for two random variables is defined as follows:66Theorem 4.3.1. For two random variables X and Y with probability distributions P and Q.Any coupling probability distribution Pˆ satisfies,||P (X = x)− P ′(Y = x)||tv ≤ 2Pˆ (Xˆ 6= Yˆ ) (4.94)For details of the proof refer [Den Hollander, 2012].4.3.1 Discrete time random walkConsider a random walk in 1 dimension X = (Xn)n≥0 taking values on state space Z startingat the origin that is X0 = 0 and Xn =∑ni=1 Yi. Here, Yi’s are i.i.d Bernoulli random variableswith,P (Yi = +1) = P (Yi = −1) = 12Theorem 4.3.2. If X is a simple random walk then for all even k ∈ Z,limn→∞||P (Xn ∈ .)− P (Xn + k ∈ .)||tv = 0Proof. LetX ′ be a independent copy ofX which starts from k, i.e., X ′0 = k and Pˆ denotesthe joint probability distribution of (X,X ′). Then, using theorem 3.2.1||P (Xn ∈ .)− P (X ′n ∈ .)||tv ≤ 2Pˆ (Xn 6= X ′n) (4.95)Next, define a random walk X˜ = (X˜n)n≥0 such that X˜n = X ′n −Xn. This random walk startswith initial position k, i.e., X˜0 = k. The random variable Y˜i denotes the increments in X˜n.Here, possible values for an increment are 1) +2X ′n goes up by one step andXn goes down byone step, 2) −2 when X ′n goes down and Xn goes up, 3) 0 when both the random walks moveup or down together. Hence distribution of Y˜i is,P˜ (Y˜i = +2) = P˜ (Y˜i = −2) = 14, P˜ (Y˜i = 0) =12. (4.96)67This is a random walk that takes steps in multiple of 2 and moves only half of the time. LetT = τ˜0 = {n ≥ 0 : X˜n = 0} (4.97), which denotes the time when Xn = X ′n. Note that the random walk X˜ is recurrent Markovchain [Durrett, 1999] which means X˜ visits state 0 infinitely often from any state. Hence recur-rence of X˜ implies that P˜ (T <∞) = 1 that means limn→∞ P˜ (Xn −X ′n = 0) = 1. Therefore,limn→∞||P (Xn ∈ .)− P (X ′n ∈ .)||tv = 0 (4.98)4.3.2 Continuous time size biased binary branching processConsider a branching process with the birth rate β and death rate µ for a particle such thatβ = µ = c. Only one particle dies at a death event, and only one offspring is born at a birthevent. Let N(t) be the population size at time t with starting population 1, i.e., N(0) = 1.Assuming that there is 0 probability of no birth or death event in a small time interval,P (N(t+∆t) = N(t) + 1|N(t)) = cN(t)∆t+ o(∆t) (4.99)P (N(t+∆t) = N(t)− 1|N(t)) = cN(t)∆t+ o(∆t) (4.100)Let N1(t) be an independent copy of process N(t), which starts with 1 particle. Now, thetransition probabilities for the size-biased version of N1(t), where bias in birth and death rateis w.r.t. the population not at the particle level.Q(N1(t+∆t) = N1(t) + 1|N1(t)) = c(N1(t) + 1)∆t+ o(∆t) (4.101)Q(N1(t+∆t) = N1(t)− 1|N1(t)) = c(N1(t)− 1)∆t+ o(∆t) (4.102)Suppose N2(t) is an independent copy of the process N(t), which starts with 2 particles, andthe size-biased version of N2(t) has similar transition probabilities as defined above.68Lemma 4.3.1. The distribution of population size for N1(t) and N2(t) as t→∞ is same,limt→∞||Q(N1(t) ∈ .)−Q(N2(t) ∈ .)||tv = 0Proof. Let Qˆ be the joint probability distribution of (N1, N2) then using theorem 3.2.1||Q(N1(t) ∈ .)−Q(N2(t) ∈ .)||tv ≤ Qˆ(N1(t) 6= N2(t)) (4.103)Next, define a process (N˜(t))t≥0 such that N˜(t) = N1(t)−N2(t). Here, N˜(0) = 1 and thepossible increments in infinitesimal time in N˜(t) are {+2,−2} with probability,Q(N˜(t+∆t) = N˜(t) + 2|N˜(t)) = c(N1(t) + 1)− c(N2(t)− 1) + o(∆t) = c(N˜(t) + 2) + o(∆t)(4.104)Q(N˜(t+∆t) = N˜(t)− 2|N˜(t)) = c(N1(t)− 1)− c(N2(t) + 1) + o(∆t) = c(N˜(t)− 2) + o(∆t)(4.105)Note that, the rate at which N˜(t) increases or decreases by 2 per unit time is identical which isc. Here, the embedded Markov chain of N˜(t) is a symmetric random walk, which has recurrentstate 0. Hence,limt→∞Q(N˜(t) = 0) = 1 =⇒ limt→∞Qˆ(N1(t)−N2(t) = 0) = 1. (4.106)Therefore Qˆ(N1(t) 6= N2(t)) = 1 and the asymptotic probability distribution of two size-biasedbranching processes with different initial state is the same.limt→∞||Q(N1(t) ∈ .)−Q(N2(t) ∈ .)||tv = 0.694.4 Reproductive value and martingaleConsider the age-structured branching process from section 3.3, where P is the probabilitymeasure of an ordinary process andQ is the probability measure for a size-biased process. Therelation (4.74) between Q and P is given by,dQxdPx∣∣∣Ft=Wx(t)Wx(0)(4.107), where Ft is the filtration of branching process till time t.Let ζ denote the extinction time, then the probability of the event thatX(t) survives till timet with respect to P measure is,∫w∈Ω1{ζ > t|Ft}dPx(w) =∫w∈ΩWx(0)1{ζ > t|Ft}dQx(w)Wx(t)(4.108)As size biased process never dies out means that 1{ζ > t} = 1 w.r.t Q measure, and we get,Px(ζ > t) =∫w∈ΩdQx(w)Wx(0)Wx(t)= EQx[Wx(0)Wx(t)], (4.109)where EQ represents expected value w.r.t Q probability measure.Theorem 4.4.1. Suppose (X1(t))t≥0 and (X2(t))t≥0 are two copies of critical age structuredbranching process (X(t))t≥0 starting from different states, i.e., x and y respectively. If repro-ductive value h is a bounded function and for all a ∈ R+ birth rate β(a) >  for some  > 0.Then there exists a coupling of X1 and X2 such that,h(x)h(y)= limt→∞W2(t)W1(t)whereWi is the martingale for the process Xi(t).Before proving Theorem 4.4.1, we need to prove two results. The first one is to show thatthere exists a successful coupling ofX1 andX2 with respect to size-biased measureQ1 andQ2.The second one is to prove the convergence of the inverse ratio of martingales w.r.t couplingmeasure Q as t→∞.70We know that in a size-biased branching process for a critical age-structured populationmodel, there exists a distinguished spinal path (section 4.2.3). Also, each particle grows itslineage in accordance with probability measure Pi if it is a non-spine particle and Qi if it is aspine particle for i ∈ {1, 2}. Now, for successful coupling ofX1 andX2, we need to couple thespine paths of both processes because if spine paths are coupled, then all the particles emanatingfrom the spine can be coupled as particles are independent of the other particles. In addition,we want lineages of non-spine particles to die out in finite time, which is growing as ordinarycritical branching processes, in which case the probability of extinction is 1. Therefore, the taskof coupling two population processes that is a set of random variables or a bunch of particlesreduces to coupling the spine paths and ensuring that non-spine branches die out.To couple the two spine paths, we need to couple the rate of spine movement. One of thespines is moving at a higher rate relative to the other one, if there is a bound on the rates suchthat the difference in rates is at most  > 0 then eventually the spines will couple.Lemma 4.4.1. For (X1(t))t≥0 and (X2(t))t≥0 age structured branching processes and ξvi(s)as the spine path for i ∈ {1, 2} in size biased version of these processes. Define the couplingtime τcpl after which both the spines move together,τcpl = inf{t : ξv1(s) = ξv2(s)}.If reproductive value h is a bounded function and birth rate β is bounded below, thenP (τcpl <∞) = 1Proof. To understand the coupling of two spine paths, let’s characterize the movement of thespine. For both the spine paths ξv1 and ξv2 , the rate at which the spine label is transferred to thenewborn is given by equation (4.76), which isλi(t) =h(0)h(ξvi(t))β(t), for i ∈ {1, 2}, (4.110)where ξvi(t) is the age of the spine particle for ith process. To couple the two spine paths, it71suffices that λi(t) > λ0 for all t ∈ R+, where λ0 > 0 because it implies that coupling time isdominated by an exponential random variable with rate λ0. We can get the lower bound λ0 bymultiplying the lower bound of the birth rate and the lower bound of the ratio h(0)/h(a). Now,lower bound on h(0)/h(a) means that for all a ∈ R+, h(0)/h(a) > h∗ for some h∗ > 0, whereh∗ =h(0)supa∈R+ h(a)(4.111), which implies that h is a bounded function and the lower bound of birth rate means that for alla ∈ R+, i.e., β(a) > β0 for some β0 > 0. Hence,λ0 = h∗β0. (4.112)This lower bound on λi(t) gives us a upper bound on coupling time τcpl, which means P (τcpl <∞) = 1. Note that, the assumptions on reproductive value function and birth rate are strong asit requires β(0) > 0, i.e., the newborn has a positive birth rate, but this assumption can be madeweaker if required.Note that, a critical age-structured branching process provides a special structure for suc-cessful coupling, i.e., the almost sure extinction of lineages of non-spine particles, which isnot possible in the case of a supercritical branching process. Hence, it is difficult to achievesuccessful coupling if the age-structured branching process is supercritical.LetWi(t) denote the martingale for the process Xi at time t given that t ≥ τcpl and decom-pose it into two parts:Wi(t) = W+i (t) +W−i (t) (4.113), where W+i denotes the total reproductive value of the spine particles of the process Xi withI ∈ {1, 2} andW−i is the total reproductive value from particles not descended from the spineat τcpl.Note that, Wx(s, t) is the normalized total reproductive potential of the process initializedwith a single particle of age x at time s, normalized to haveE[Wx(s, t)] = 1 for all t > 0, which72is achieved by takingWx(s, t) =Wx(t)Wx(s)=Wx(t)h(x), (4.114)so that Wx(s, s) = 1. In the above expression, second equality holds as at time s, there isonly one particle of age x, i.e., Wx(s) = h(x). Here, W follows the semigroup property, i.e.,Wx(s, t+ k) = Wx(s, t)Wx(t, t+ k) as shown below,Wx(s, t+ k) =Wx(t+ k)Wx(s)(4.115)=Wx(t+ k)Wx(s)Wx(t)Wx(t)=Wx(t)Wx(s)Wx(t+ k)Wx(t)(4.116)= Wx(s, t)Wx(t, t+ k). (4.117)Let’s compute the instantaneous drift of 1/W w.r.t. Q measure as a function of the popu-lation measure and let’s denote it by νQ(X), where X is the finite counting measure over theset of type/state space. The drift of 1/W consists of the deterministic part due to aging andthe random jump part. The former is unaffected by size biasing, while the latter is affectedby multiplying rates by the relative values of the martingale being used for size-biasing. WithX as the counting measure, W = h(X) =∑x∈X h(x). Since age increases at deterministicspeed 1 while alive, using the chain rule, the deterministic part of the drift of 1/W = 1/h(X)is −∑x∈X h′(x)/h2(X). Let ∆ be counting measure which denotes the change in value ofWat time t, total drift of 1/W is,νQ(X) =−1h2(X)∑x∈Xh′(x) +∑∆( 1h(X) + ∆− 1h(X))h(X) + ∆h(X)Q(X,X +∆) (4.118)=−1h2(X)∑x∈Xh′(x) +∑∆−∆h2(X)Q(X,X +∆) (4.119)=−1h2(X)(∑x∈Xh′(x) +∑∆∆Q(X,X +∆)), (4.120)where∑x∈X means summation over finitely many particle’s states and Q(X,X + ∆) is therate at which the counting measure X changes by ∆. Note that, νP (W ) =∑x∈X h′(x) +73Figure 4.1: Coupling of two size-biased age-structured branching processes starting with agesx (blue path) and y (green path). At the time τcpl, the spine path of both the processes coupleafter which the red path represents the trajectory. After coupling, the blue and green branchesrepresent lineages of uncoupled particles, and the black cross represents the termination of un-coupled lineages. The contribution of blue and green branches after coupling is W−i , and thecontribution of red branches isW+i .∑∆∆Q(X,X + ∆), i.e, the drift of W w.r.t P measure, which is equal to 0 since W is a P -martingale. From (4.120), we can conclude that the drift of 1/W is negativewhen h(X)+∆ = 0and the drift is 0when h(X)+∆ 6= 0, with respect toQmeasure. Hence, 1/W is a non-negativeQ-supermartingale, then usingmartingale convergence theorem [Durrett, 2019], we can say that1/W converges almost surely, i.e.,limt→∞1W (t)=1W∞(4.121)such that E[1/W∞] <∞.Lemma 4.4.2. For (X1(t))t≥0 and (X2(t))t≥0 critical age-structured branching processes withξvi(s) as the spine path for i ∈ {1, 2} if P (τcpl <∞) = 1 then,limt→∞W2(t)W1(t)= limt→∞W+2 (t)W+1 (t), (4.122)74where for t > τcpl,W+i (t) is given by,W+i (t) =h(xτcpl)h(xi)Wxcpl(τcpl, t). (4.123)Here, xτcpl is the state of spine particle at coupling time τcpl.Proof. Let’s estimate the upper bound for the ratio of the martingales for t ≥ τcpl,W2(t)W1(t)=W+2 (t) +W−2 (t)W+1 (t) +W−1 (t)≤ W+2 (t) +W−2 (t)W+1 (t)=(W+2 (t)W+1 (t)+W−2 (t)W+1 (t))(4.124)Note that, non-spine particles for X1 and X2 follow the law of ordinary critical branch-ing processes, so their lineages go extinct almost surely for large t w.r.t P measure, i.e.,limt→∞ P (ζ1 > t) = 0. Thus the total reproductive value of non-spine particle of X2 pro-cess,W−2 (t)→ 0 as t→∞ and using the fact that limt→∞ 1/W+1 = 1/W∞ <∞, implies thatlimt→∞W−2 (t)/W+1 (t) = 0. Assuming that limt→∞W+2 (t)/W+1 (t) exists,limt→∞supW2(t)W1(t)≤ limt→∞W+2 (t)W+1 (t). (4.125)Next, estimate the lower bound for the ratio of martingales,W2(t)W1(t)=W+2 (t) +W−2 (t)W+1 (t) +W−1 (t)≥ W+2 (t)W+1 (t) +W−1 (t)≥ W+2 (t)W+1 (t)(1− W−1 (t)W+1 (t))(4.126)Again, using the fact that W−1 (t) → 0 as t → ∞ and limt→∞ 1/W+2 = 1/W∞ < ∞, impliesthat limt→∞W−1 (t)/W+1 (t) = 0. Assuming that limt→∞W+2 (t)/W+1 (t) exists,limt→∞infW2(t)W1(t)≥ limt→∞W+2 (t)W+1 (t)(4.127)This concludes the proof of lemma 4.4.2. Next is the proof of theorem 4.4.1.Proof. We know that W+i represents the total reproductive value of the particles descendentfrom coupled spine paths after and at coupling time. Now, for t > τcpl, W+i (t) is productof reproductive value at coupling time (which is same for both processes) relative to initial75reproductive value of ith process and total reproductive value of the coupled spine particlesfrom τcpl till time t.W+i (t) =h(xτcpl)h(xi)Wxcpl(τcpl, t), (4.128)where xcpl is the state of spine particle at τcpl which is same for both the spine paths, i.e., xcpl =ξv1(τcpl) = ξv2(τcpl). Hence, using lemma 4.4.2 and (4.128),limt→∞W2(t)W1(t)= limt→∞W+2 (t)W+1 (t)=h(x1)h(x2)=h(x)h(y)(4.129)76Chapter 5ConclusionIn this work, we utilized two established methods, namely PF theory and renewal theory, to de-fine the reproductive value function for various deterministic and stochastic population models.We derived two new characterizations of reproductive value. The first main result demonstratesthat the limiting ratio of survival probability is equal to the ratio of initial reproductive valuesfor two critical age-structured population models that commence from different state-time loca-tions. The second result establishes the existence of a successful coupling between size-biasedversions of two critical age-structured population models, leading to the equality of the ratio ofinitial reproductive values with the limit of the inverse ratio of martingales. The new charac-terizations of reproductive value can be utilized in population genetics to compare the survivalprobability of genes by estimating the reproductive value of their hosts.Note that, in the first characterization, we estimate the rate at which survival probabilitydecays as the process goes extinct almost surely, whereas, in the second characterization ofreproductive value, we are comparing two age-structured branching processes while they arealive. For the second characterization, the martingale property of reproductive potential is ofimportance, sincewe need amartingale to construct a size-biased version of the ordinary branch-ing process, which is required for the successful coupling of the two processes to compare theprocesses while they are alive. For the successful coupling of two branching processes, we haveto couple a collection of particles from the two processes, which is a challenging task becausethe particles follow random processes for reproduction and survival, and size-biasing reduces77the task of coupling the collection of particles to coupling the unique spine path.Our analysis focused on single-type and multi-type population models with state-dependentbirth and death rates parameters. In the context of deterministic models, we defined reproduc-tive value for linear population models. Additionally, we demonstrated the fundamental prop-erty of reproductive value as an invariant function for time-homogeneous population models. Inthe case of stochastic population models, we established the concept of reproductive potentialfor both single and multi-type population models. Using the well-known criteria of dominanteigenvalue, we classify the stochastic population processes into three categories: sub-critical,critical, and super-critical process. Notably, we demonstrated that the reproductive potentialexhibits the property of being a martingale in critical stochastic population models, which sig-nifies its time invariance.In the study of population models, branching process formulation provides a rich form ofthe process as compared to their deterministic counterparts. Specifically, random demographicand environmental fluctuations have a significant impact on extinction events, which is difficultto capture in deterministic models. In the case of processes that start with a few individuals,stochasticity plays an important role in deriving the evolution of the process. For example, atthe start of an epidemic, it is useful to study the individual-level dynamics to understand thetrajectory of the spread.In epidemic modeling, our work can be applied to understand the disease dynamics, whereit is reasonable to assume that the transmission rate and recovery rate are the same within anage group. We can estimate the reproductive value, and use it to get the extinction probabilityof the disease w.r.t the initial state. In the study of cell cycles, where at the end of a life cycle,the cell doubles in size and divides into two daughter cells of approximately the same size, theprocess can be formulated as a branching process with the independence assumption. Using thereproductive value characterizations, the long-term behavior of these types of processes can beanalyzed and survival probability can be estimated. In population genetics, genetic mutationsurvival is important to access the probability of fixation or extinction in the long run. Our workcan be applied to estimate the survival probability using the reproductive value characterization.In addition, using themartingale change ofmeasure approach, we can compare the gene survival78w.r.t different initial environmental conditions. An interesting example of this type of branchingprocess is the HIV mutation within an individual [Olabode et al., 2021].In future work, it may be possible to use the second characterization to define the repro-ductive value in time inhomogeneous branching population models where birth and death ratesdepend on both time and the state of the particle. In such cases, the application of renewal theorybecomes challenging because, despite the birth age being zero for all particles, the lifespan andreproduction process differ for particles born at different times. Consequently, the distributionsof lifespan and reproduction process no longer exhibit identical distributions.Specifically, for the time inhomogeneous case, we are trying to obtain a function similarto the reproductive value function, i.e., a martingale. To define this martingale, we need tointroduce compensators, which are very challenging to deal with. These compensators representthe inverse of the expected growth rate along the trajectory of an individual. There are two keychallenges associated with the compensator. Firstly, the compensator tends to have small valueswhen an individual reproduces later in its lifetime. This issue can be addressed by imposingan upper bound on the birth rate, thereby controlling the compensator’s behavior. Secondly,the dependence of the compensator on the initial state of the branching process complicates theanalysis.Further, in the future, it may be possible to extend our characterizations for the populationprocess in which individuals interact with others, i.e., the reproduction and survival processesare dependent on the other individuals in the population. One approach for formulating theseprocesses is explored in [Diekmann et al., 2003] for steady-state analysis, where individuals areconsidered independent, and the dependency conditions are taken as the environmental condi-tions with a consistent feedback loop connecting the individuals with the environment.79Bibliography[Al-Arydah et al., 2012] Al-Arydah, M., Lutscher, F., et al. (2012). Modeling gender-structured wildlife diseases with harvesting: Chronic wasting disease as an example. In-ternational Scholarly Research Notices, 2012.[Arratia et al., 2019] Arratia, R., Goldstein, L., and Kochman, F. (2019). Size bias for one andall. Project Euclid.[Athreya and Ney, 1982] Athreya, K. B. and Ney, P. (1982). A renewal approach to the perron-frobenius theory of non-negative kernels on general state spaces. Mathematische Zeitschrift,179(4):507–529.[Athreya and Ney, 2004] Athreya, K. B. and Ney, P. E. (2004). Branching processes. CourierCorporation.[Auger et al., 2008] Auger, P., Magal, P., and Ruan, S. (2008). Structured population modelsin biology and epidemiology, volume 1936. Springer.[Belomestny et al., 2022] Belomestny, D., Moulines, E., and Samsonov, S. (2022). Variancereduction for additive functionals of markov chains via martingale representations. Statisticsand Computing, 32(1):16.[Berryman, 1992] Berryman, A. A. (1992). The orgins and evolution of predator-prey theory.Ecology, 73(5):1530–1535.[Bianca and Dogbe, 2017] Bianca, C. and Dogbe, C. (2017). On the existence and unique-ness of invariant measure for multidimensional stochastic processes. Nonlinear Studies-TheInternational Journal.80[Biggins and Kyprianou, 2004] Biggins, J. D. and Kyprianou, A. E. (2004). Measure changein multitype branching. Advances in Applied Probability, 36(2):544–581.[Boenkost et al., 2022] Boenkost, F., Foutel-Rodier, F., and Schertzer, E. (2022). The ge-nealogy of a nearly critical branching processes in varying environment. arXiv preprintarXiv:2207.11612.[Brodie and Post, 2010] Brodie, J. F. and Post, E. (2010). Nonlinear responses of wolverinepopulations to declining winter snowpack. Population Ecology, 52:279–287.[Cloez et al., 2021] Cloez, B., de Saporta, B., and Roget, T. (2021). Long-time behavior anddarwinian optimality for an asymmetric size-structured branching process. Journal of Math-ematical Biology, 83:1–30.[Clutton-Brock et al., 1997] Clutton-Brock, T., Rose, K., and Guinness, F. (1997). Density–related changes in sexual selection in red deer. Proceedings of the Royal Society of London.Series B: Biological Sciences, 264(1387):1509–1516.[Cox et al., 2022] Cox, A. M., Horton, E., and Villemonais, D. (2022). Binary branching pro-cesses with moran type interactions. arXiv preprint arXiv:2207.03323.[Cui and Chen, 2007] Cui, M. and Chen, Z. (2007). The exact solution of nonlin-ear age-structured population model. Nonlinear Analysis: Real World Applications,8(4):1096–1112.[Cushing, 1988] Cushing, J. (1988). Nonlinear matrix models and population dynamics. Nat-ural Resource Modeling, 2(4):539–580.[Den Hollander, 2012] Den Hollander, F. (2012). Probability theory: The couplingmethod. Lecture notes available online (http://websites. math. leidenuniv. nl/probability/lec-turenotes/CouplingLectures. pdf).[Diekmann et al., 2003] Diekmann, O., Gyllenberg, M., andMetz, J. (2003). Steady-state anal-ysis of structured population models. Theoretical population biology, 63(4):309–338.81[Doebeli and Koella, 1994] Doebeli, M. and Koella, J. C. (1994). Sex and population dy-namics. Proceedings of the Royal Society of London. Series B: Biological Sciences,257(1348):17–23.[Durrett, 1999] Durrett, R. (1999). Essentials of stochastic processes, volume 1. Springer.[Durrett, 2019] Durrett, R. (2019). Probability: theory and examples, volume 49. Cambridgeuniversity press.[Feller, 1941] Feller, W. (1941). On the integral equation of renewal theory. The Annals ofMathematical Statistics, 12(3):243–267.[Feller, 1991] Feller, W. (1991). An introduction to probability theory and its applications,Volume 2, volume 81. John Wiley & Sons.[Fisher, 1999] Fisher, R. A. (1999). The genetical theory of natural selection: a completevariorum edition. Oxford University Press.[Foucart and Vidmar, 2022] Foucart, C. and Vidmar, M. (2022). Continuous-state branchingprocesses with collisions: first passage times and duality. arXiv preprint arXiv:2207.04621.[Foutel-Rodier and Schertzer, 2022] Foutel-Rodier, F. and Schertzer, E. (2022). Convergenceof genealogies through spinal decomposition, with an application to population genetics.arXiv preprint arXiv:2201.12412.[Frobenius et al., 1912] Frobenius, G., Frobenius, F. G., Frobenius, F. G., Frobenius, F. G.,and Mathematician, G. (1912). Uber matrizen aus nicht negativen elementen. KoniglicheAkademie der Wissenschaften Sitzungsber, Kon.[Fujiwara and Caswell, 2001] Fujiwara, M. and Caswell, H. (2001). Demography of the en-dangered north atlantic right whale. Nature, 414(6863):537–541.[Green, 1977] Green, P. (1977). Conditional limit theorems for general branching processes.Journal of Applied Probability, 14(3):451–463.82[Guzowska et al., 2011] Guzowska, M., Luis, R., and Elaydi, S. (2011). Bifurcation and in-variant manifolds of the logistic competition model. Journal of Difference Equations andApplications, 17(12):1851–1872.[Hardy and Harris, 2006] Hardy, R. and Harris, S. C. (2006). Spine proofs for lp-convergenceof branching-diffusion martingales. arXiv preprint math/0611056.[Heathcote et al., 1967] Heathcote, C. R., Seneta, E., and Vere-Jones, D. (1967). A refine-ment of two theorems in the theory of branching processes. Theory of Probability and ItsApplications, 12(2):297–301.[Henson and Cushing, 1997] Henson, S. M. and Cushing, J. (1997). The effect of periodichabitat fluctuations on a nonlinear insect population model. Journal of Mathematical Biol-ogy, 36(2):201–226.[Hong et al., 2022] Hong, W., Liang, S., and Zhang, X. (2022). Conditional-convergence forthe martingale of a critical branching process in random environment. Proceedings of theSteklov Institute of Mathematics, 316(1):184–194.[Jagers, 1989] Jagers, P. (1989). General branching processes as markov fields. StochasticProcesses and their Applications, 32(2):183–212.[Jagers and Klebaner, 2000] Jagers, P. and Klebaner, F. C. (2000). Population-size-dependentand age-dependent branching processes. Stochastic Processes and their Applications,87(2):235–254.[Joffe and Spitzer, 1967] Joffe, A. and Spitzer, F. (1967). On multitype branching processeswith �� 1. Journal of mathematical analysis and applications, 19(3):409–430.[Jorgensen and Gromiec, 2016] Jorgensen, S. and Gromiec, M. (2016). Developments in envi-ronmental modelling. Elsevier.[Kesten et al., 1966] Kesten, H., Ney, P., and Spitzer, F. (1966). The galton-watson pro-cess with mean one and finite variance. Theory of Probability and Its Applications,11(4):513–540.83[Kesten and Stigum, 1966] Kesten, H. and Stigum, B. P. (1966). A limit theorem formultidimensional galton-watson processes. The Annals of Mathematical Statistics,37(5):1211–1223.[Knowlton, 2004] Knowlton, N. (2004). Multiple “stable” states and the conservation of ma-rine ecosystems. Progress in Oceanography, 60(2-4):387–396.[Kolesko and Sava-Huss, 2023] Kolesko, K. and Sava-Huss, E. (2023). Limit theorems fordiscrete multitype branching processes counted with a characteristic. Stochastic Processesand their Applications, 162:49–75.[Krein and Rutman, 1948] Krein, M. G. and Rutman, M. A. (1948). Linear operators leavinginvariant a cone in a banach space. Uspekhi mat. nauk, 3(1):3–95.[Kurtz et al., 1997] Kurtz, T., Lyons, R., Pemantle, R., and Peres, Y. (1997). A conceptualproof of the kesten-stigum theorem for multi-type branching processes. Classical and mod-ern branching processes, pages 181–185.[Kyprianou, 2000] Kyprianou, A. E. (2000). Martingale convergence and the stopped branch-ing random walk. Probability theory and related fields, 116(3):405–419.[Kyprianou and Palmowski, 2005] Kyprianou, A. E. and Palmowski, Z. (2005). A martingalereview of some fluctuation theory for spectrally negative lévy processes. Séminaire de Prob-abilités XXXVIII, pages 16–29.[Lande et al., 2006] Lande, R., Engen, S., Sther, B.-E., and Coulson, T. (2006). Estimatingdensity dependence from time series of population age structure. The American Naturalist,168(1):76–87.[Leslie, 1945] Leslie, P. H. (1945). On the use of matrices in certain population mathematics.Biometrika, 33(3):183–212.[Li and Schneider, 2002] Li, C.-K. and Schneider, H. (2002). Applications of perron–frobeniustheory to population dynamics. Journal of mathematical biology, 44(5):450–462.84[Li et al., 2019] Li, P.-S., Yang, X., and Zhou, X. (2019). A general continuous-state nonlinearbranching process.[Liggett et al., 1999] Liggett, T. M. et al. (1999). Stochastic interacting systems: contact, voterand exclusion processes, volume 324. springer science & Business Media.[Lion, 2017] Lion, S. (2017). On the dynamics of reproductive values and phenotypic traits inclass-structured populations. bioRxiv, page 155879.[Longstaff, 1977] Longstaff, B. C. (1977). The dynamics of collembolan populations: a matrixmodel of single species population growth. Canadian Journal of Zoology, 55(2):314–324.[Lotka, 1939a] Lotka, A. J. (1939a). A contribution to the theory of self-renewing aggregates,with special reference to industrial replacement. The Annals of Mathematical Statistics,10(1):1–25.[Lotka, 1939b] Lotka, A. J. (1939b). On an integral equation in population analysis. The Annalsof Mathematical Statistics, 10(2):144–161.[Lyons et al., 1995] Lyons, R., Pemantle, R., and Peres, Y. (1995). Conceptual proofs of llog l criteria for mean behavior of branching processes. The Annals of Probability, pages1125–1138.[Malthus and Winch, 1992] Malthus, T. R. and Winch, D. (1992). An essay on the principle ofpopulation. Cambridge University Press Cambridge.[Manly, 1990] Manly, B. F. (1990). Stage-structured populations. Springer.[M’kendrick, 1925] M’kendrick, A. (1925). Applications of mathematics to medical problems.Proceedings of the Edinburgh Mathematical Society, 44:98–130.[Moharrami et al., 2020] Moharrami, M., Subramanian, V., Liu, M., and Sundaresan,R. (2020). The erlang weighted tree, a new branching process. arXiv preprintarXiv:2002.03993.85[Molchanov and Yarovaya, 2012] Molchanov, S. and Yarovaya, E. (2012). Branching pro-cesses with lattice spatial dynamics and a finite set of particle generation centers. InDokladyMathematics, volume 86, pages 638–641. SP MAIK Nauka/Interperiodica.[Murray, 2002] Murray, J. D. (2002). Mathematical biology: I. An introduction. Springer.[Ojeda and Pardo, 2020] Ojeda, G. B. and Pardo, J. C. (2020). Branching processes with pair-wise interactions. arXiv preprint arXiv:2009.11820.[Olabode et al., 2021] Olabode, D., Rong, L., and Wang, X. (2021). Stochastic investigationof hiv infection and the emergence of drug resistance. Mathematical Biosciences and Engi-neering, 19(2).[Pavliotis and Zanoni, 2022] Pavliotis, G. A. and Zanoni, A. (2022). Eigenfunction martin-gale estimators for interacting particle systems and their mean field limit. SIAM Journal onApplied Dynamical Systems, 21(4):2338–2370.[Perron, 1907] Perron, O. (1907). Neue kriterien fur die irreduzibilitat algebraischer gleichun-gen. Walter de Gruyter, Berlin/New York Berlin, New York.[Rose et al., 2005] Rose, K. E., Louda, S. M., and Rees, M. (2005). Demographic and evo-lutionary impacts of native and invasive insect herbivores on cirsium canescens. Ecology,86(2):453–465.[Ruohomaki et al., 2003] Ruohomaki, K., Klemola, T., Kaitaniemi, P., and Kaar, M. (2003).Crowding-induced responses in a geometrid moth revisited: a field experiment. Oikos,103(3):489–496.[Silva Matos et al., 1999] Silva Matos, D. M., Freckleton, R. P., and Watkinson, A. R. (1999).The role of density dependence in the population dynamics of a tropical palm. Ecology,80(8):2635–2650.[Smith et al., 1977] Smith, D. P., Keyfitz, N., and Lotka, A. J. (1977). The stability of thenormal age distribution. Mathematical Demography: Selected Papers, pages 101–107.86[Solbrig et al., 1988] Solbrig, O. T., Sarandon, R., and Bossert, W. (1988). A density-dependent growth model of a perennial herb, viola fimbriatula. The American Naturalist,131(3):385–400.[Tanner, 1999] Tanner, C. (1999). Constraints on environmental behaviour. Journal of envi-ronmental psychology, 19(2):145–157.87
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Foxall, Eric
Reproductive value for branching population models
Text