Kinetic theory for structured populations: application to stochastic sizer-timer models of cell proliferation
KKinetic theory for structured populations: application to stochastic sizer-timermodels of cell proliferation
Mingtao Xia and Tom Chou
1, 2 Department of Mathematics, UCLA, Los Angeles, CA, 90095-1555, USA Department of Computational Medicine, UCLA, Los Angeles, CA, 90095-1766, USA ∗ We derive the full kinetic equations describing the evolution of the probability density distributionfor a structured population such as cells distributed according to their ages and sizes. The kineticequations for such a “sizer-timer” model incorporates both demographic and individual cell growthrate stochasticities. Averages taken over the densities obeying the kinetic equations can be usedto generate a second order PDE that incorporates the growth rate stochasticity. On the otherhand, marginalizing over the densities yields a modified birth-death process that shows how ageand size influence demographic stochasticity. Our kinetic framework is thus a more complete modelthat subsumes both the deterministic PDE and birth-death master equation representations forstructured populations.
Keywords : Age Structure, Birth-Death Process, Kinetics, Fission
I. INTRODUCTION
Across many diverse applications, mathematical models have been formulated to describe the evolution of popula-tions according to a number of individual attributes such as age, size, and/or added size since birth. For example,deterministic age-structured models that incorporate age-dependent birth and death were developed by McKendrickand have been applied to human populations [1]. More recently, there has been renewed interest in cell size control[2, 3], cellular division mechanisms [4], and structured cell population models [5, 6].When considering proliferating cell populations, individual cell growth is interrupted by cell division events thatgenerate smaller daughter cells. Cell division is a process that involves many biochemical steps and complex biophysicalmechanisms that involves metabolism, gene expression, protein production, DNA replication, chromosome separation(for eukaryotic cells), and fission or cell wall formation [7–11]. To simplify the understanding of which factors triggercell division, three basic models that subsume these complex processes have been proposed. Cells can divide basedon their age since birth, volume (size), or added volume since birth y [2, 12]. PDE approaches for the timer, sizer,and adder models, as well as combinations of these models, have been well-studied [6, 13, 14]. These PDE approachesimplicitly describe the mean density of cells in age, size, and/or added size, and are considered deterministic models.However, there has been much less development of structured populations models that incorporate stochastic effects.In the presence of stochasticity, how would the PDEs be modified? In the sizer-timer type of structured populationmodels, stochasticity can arise in the growth dynamics of each cell as well as in the random times of cell division anddeath (demographic stochasticity).Stochasticity arising from random times of birth and death (demographic stochasticity) has been considered in timer-like models for age-structured populations [15, 16]. This approach generalized the classic deterministic McKendrickequation to a higher dimension (dynamically varying) associated with the number of individuals in the system. Thishigher-dimensional stochastic “kinetic theory” allows one to systematically connect an age-indepedent birth-deathmaster equation description to the deterministic age-structured McKendrick model. A comprehensive and generaltreatment of the age-structured stochastic process using a Doi-Peliti operator formulism has also been developed forcalculation of correlation functions [17]. The full kinetic theory has only been developed for age-structured populationsand only includes demographic stochasticity (since chronological age is a deterministic quantity proportional to time).Other approaches using stochastic hybrid systems [18] have been used to incorporate the influence of random birthtimes of population-level variations in cell size. Intrinsic stochasticity in the growth rate of an individual cell has beentreated in terms of Langevin equations for cell size [19], effective potentials [3] and stochastic maps [12, 20]. Recently,Chapman-Kolmogorov equations have also been applied to study the effect of different sources of noise in cellularproliferation [21]. However, stochasticity in the intrinsic growth rate has not been considered within demographicallystochastic kinetic theory. ∗ Electronic address: [email protected] a r X i v : . [ q - b i o . P E ] J a n In this paper, we shall derive a kinetic theory for the sizer-timer model of cell proliferation that incorporates bothdemographic stochasticity and intrinsic stochasticity in the growth of individual cells. In the next section, we derive theFokker-Planck equation for the size of an individual cell and define the probabilistic quantities needed to construct thefull kinetic theory. This equation is then marginalized in Section III to explicitly isolate and show the feature limits ofintrinsic stochasticity and demographic stochasticity. Including both sources of stochasticity renders the calculationsof marginalized densities rather technical, but by successively taking the marginalized single-density limits, we showhow the theory reduces to simpler forms and reveal the procedure for solving the full high-dimensional problem.Moreover, by taking higher moments of the density, an unclosed hierarchy of equations that reflect demographicstochasticity arises. Our results generalize a large body of work on sizer-timer PDE models to include stochasticprocesses, both at the individual and population levels.
II. DERIVATION OF KINETIC THEORY
Here we outline the derivation of the kinetic equation for a population of dividing cells of different ages a and sizes(volumes) x . We start from the SDE for the size of a single cell at time t :d X t = g ( X t , A t , t )d t + σ ( X t , A t , t )d W t , X t , A t ∈ Λ , (1)where Λ := [0 , ∞ ), A t is the cell’s age (time that has elapsed after its birth), g ( X t , A t , t ) is the size- and age-dependentgrowth rate, and W t is a standard Wiener process with independent, normally distributed increments W t − W s , zeromean, and variance t − s . The parameter σ ( X t , A t , t ) represents the strength of stochasticity in cell’s growth rate.Here, we assume both g and σ are Lipschitz continuous to ensure the existence and uniqueness of X t given any initialconditions X > , A ≥
0. We also assume σ ∈ C , σ (0 , t, a ) = ∂ x σ (0 , t, a ) = 0 so that the noise vanishes at x = 0and X t remains positive.Next, we investigate a system of m + 2 n cells, where m is the number of individual cells (singlets) and n is thenumber of twins (doublets). A twin means two daughter cells generated from the division of a common mother cell,and therefore they have the identical age. In this section, we use the notation X ( m ) t = ( X t , X t , ..., X mt ) , Y (2 n ) t = ( Y t , ..., Y nt ) , A ( m ) t = ( A t , A t , ..., A mt ) , B ( n ) t = ( B t , ..., B nt ) , (2)where A ( m ) t and B ( n ) t are ordered ages such that A it ≥ A jt ≥ , B it ≥ B jt ≥ , ∀ i > j and X ( m ) t and Y (2 n ) t are thevectors of the volumes of the m singlets and 2 n doublets that are of ages A ( m ) t and B ( n ) t , respectively, at time t . Notethat two cells in a doublet have the same age but can have different sizes; thus, the age vector B ( n ) t of the 2 n twinsstores n ages, while the size vector Y (2 n ) t stores 2 n sizes.Formally solving Eq. (1), each X it and Y jt satisfies X it = X it (cid:48) + (cid:90) tt (cid:48) g ( X is , A is , s )d s + (cid:90) tt (cid:48) σ ( X s , A s , s )d W is ,Y jt = Y jt (cid:48) + (cid:90) tt (cid:48) g ( Y js , B [ j +12 ] s , s )d s + (cid:90) tt (cid:48) σ ( Y js , B [ j +12 ] s , s )d W m + js , (3)where d W is , d W m + js are the intrinsic, indepedent fluctuations in growth rates. We assume that cell division rates areregulated by a “timer” mechanism and does not depend on cell size, i.e. , the probability that a cell in a populationof m singlets and n doublets divides during ( t, t + ∆ t ] is β m,n ( A t , t )d t + o (d t ), a function of its age A t , time t andpopulation sizes m, n . The mathematical analyses that follow require that the birth rate is independent of a cell’ssize X t . Finally, we take the continuous time limit and assume that in a finite number of cells, the possibility of twocells dividing in ( t, t + d t ] is o (d t ) as d t → Alternatively, X t might also represent the log of the cell size A. The forward equation
We evaluate the increment in time by Ito’s formula applied to a function f m,n ( X ( m ) t , Y (2 n ) t , t ; A ( m ) t (cid:48) , B ( n ) t (cid:48) ) of m in-dividual and n twin sizes given initial sizes and ages A ( m ) t (cid:48) , B ( n ) t (cid:48) at t (cid:48) < t , where the ages are defined to be in thedescending order A ≥ A ... ≥ A m ≥ B ≥ B ... ≥ B n ≥
0. Ordering the ages allows us to easily incorporate celldivision as a boundary condition in which newborn cells are represented by B n = 0. f m,n ( X ( m ) t +d t , Y (2 n ) t +d t , t + d t ; A ( m ) t (cid:48) , B ( n ) t (cid:48) ) − f m,n ( X ( m ) t , Y (2 n ) t , t ; A ( m ) t (cid:48) , B ( n ) t (cid:48) ) (cid:90) t +d tt ∂f m,n ∂s + m (cid:88) i =1 g ( X is , A is , s ) ∂f m,n ∂X is + n (cid:88) j =1 g ( Y js , B [( j +1) / s , s ) ∂f m,n ∂Y js + 12 m (cid:88) i =1 σ ( X is , A is , s ) ∂ f m,n ( ∂X is ) + 12 n (cid:88) j =1 σ ( Y js , B [( j +1) / s , s ) ∂ f m,n ( ∂Y js ) d s + m (cid:88) i =1 (cid:90) t +d tt σ ( X is , A is , s ) ∂f m,n ∂X is d W is + n (cid:88) j =1 (cid:90) t +d tt σ ( Y js , B [( j +1) / s , s ) ∂f m,n ∂Y js d ˜ W js . (4)After taking the expectation of Eq. (4) we find E [ f m,n ( X ( m ) t +d t , Y (2 n ) t +d t , t + d t ; A ( m ) t (cid:48) , B ( n ) t (cid:48) )] − E [ f m,n ( X ( m ) t , Y (2 n ) t , t ; A ( m ) t (cid:48) , B ( n ) t (cid:48) )] = E (cid:90) t +d tt d s ∂f m,n ∂s + m (cid:88) i =1 g ( X is , A is , s ) ∂f m,n ∂X is + n (cid:88) j =1 g ( Y js , B [( j +1) / s , s ) ∂f m,n ∂Y js + 12 m (cid:88) i =1 ∂ f m,n ( ∂X is ) σ ( X is , A is , s ) + 12 n (cid:88) j =1 ∂ f m,n ( ∂Y js ) σ ( Y js , B [( j +1) / s , s ) . (5)Specifically, we can take f m,n in Eq. (5) as a distribution of the form f m,n ( X ( m ) t , Y (2 n ) t , t ; A ( m ) t (cid:48) , B ( n ) t (cid:48) ) = m (cid:89) i =1 δ ( X i − X it ) n (cid:89) j =1 δ ( Y j − Y jt ) S ,m ( t ; t (cid:48) , A ( m ) t (cid:48) ) S ,n ( t ; t (cid:48) , B ( m ) t (cid:48) ) , (6)where S ,m and S ,n are joint survival possibilities S ,m ( t ; t (cid:48) , A ( m ) ) = m (cid:89) i =1 e − (cid:82) tt (cid:48) β m,n ( A i − t (cid:48) + s,s )ds , S ,n ( t ; t (cid:48) , B ( n ) ) = n (cid:89) j =1 ( e − (cid:82) tt (cid:48) β m,n ( B j − t (cid:48) + s,s )ds ) , (7)where the birth rate β ≡ β m,n can implicitly depend on the populations m, n .Next, we define ˆ p ( X ( m ) t , Y (2 n ) t , t | X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) ) as the probability density of m singlets of volumes X ( m ) t and n doublets of volumes Y (2 n ) t at time t , conditioned on there being m singlets of volumes X ( m ) t (cid:48) and ages A ( m ) t (cid:48) and n doublets with volumes Y (2 n ) t (cid:48) and ages B (2 n ) t (cid:48) at time t (cid:48) , and that no cell division occurs during [ t (cid:48) , t ]. Thequantity ˆ p ( X ( m ) t , Y (2 n ) t , t | X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) ) S ,m ( t ; t (cid:48) , A ( m ) t (cid:48) ) S ,n ( t ; t (cid:48) , B ( m ) t (cid:48) ) is thus the probability measure that thecell population at time t contains m singlets of size X ( m ) t and n doublets of size Y ( n ) t with no cell division occurringwithin [ t (cid:48) , t ], conditioned on it containing m singlets with volumes X ( m ) t (cid:48) and ages A ( m ) t (cid:48) and n doublets with volumes Y (2 n ) t (cid:48) and ages B ( n ) t (cid:48) at t (cid:48) .After substitution of the f m,n defined in Eq. (6) into Eq. (5), dividing by d t , and taking the d t → ∂∂t (cid:16) ˆ p ( X ( m ) , Y (2 n ) , t | X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) ) S ,m ( t ; t (cid:48) , A ( m ) t ) S ,n ( t ; t (cid:48) , B ( n ) t ) (cid:17) = (cid:90) Λ m d X ( m ) t (cid:90) Λ n d Y (2 n ) t ˆ p ( X ( m ) t , Y (2 n ) t , t | X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) ) (cid:20) ∂f∂t + m (cid:88) i =1 g ( X it , A it , t ) ∂f∂X it + n (cid:88) j =1 g ( Y jt , B [( j +1) / t , t ) ∂f∂Y jt + 12 m (cid:88) i =1 ∂ f∂ ( X it ) σ ( X it , A it , t ) + 12 n (cid:88) j =1 ∂ f∂ ( Y jt ) σ ( Y jt , B [( j +1) / t , t ) (cid:21) = − (cid:34)(cid:18) m (cid:88) i =1 β m,n ( A it , t ) + 2 n (cid:88) j =1 β m,n ( B jt , t ) (cid:19) ˆ p m,n + m (cid:88) i =1 ∂ (ˆ pg ( X it , A it , t )) ∂X it + n (cid:88) j =1 ∂ (ˆ pg ( Y jt , B [( j +1) / t , t )) ∂Y jt − m (cid:88) i =1 ∂ (ˆ pσ ( X it , A it , t ))( ∂X it ) − n (cid:88) j =1 ∂ (ˆ pσ ( Y jt , B jt , t ))( ∂Y jt ) (cid:35) S ,m S ,n , (8)where the last equality arises from integration by parts.Finally, we derive the PDE satisfied by the unconditioned probability density p m,n ( X ( m ) t , Y (2 n ) t , A ( m ) t , B ( n ) t , t ) given p m,n ( X ( m ) , Y (2 n ) , A ( m ) , B ( n ) , t (cid:48) ). First, we note that if no division has occurred in [ t (cid:48) , t ] and t − t (cid:48) < min { A ( m ) t , B ( n ) t } ,a system at t with m singlets of volumes X ( m ) t and ages A ( m ) t and n doublets with volumes Y (2 n ) t and ages B ( n ) t can result only from a system at t (cid:48) with m singlets with ages A ( m ) t (cid:48) = A ( m ) t − ( t − t (cid:48) ) and n doubletswith ages B ( n ) t (cid:48) = B ( n ) t − ( t − t (cid:48) ). Thus, we use the Chapman-Kolmogorov relation between the two quantitiesˆ p ( X ( m ) t , Y (2 n ) t , t | X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) ) S ,m ( t ; t (cid:48) , A ( m ) t (cid:48) ) S ,n ( t ; t (cid:48) , B ( m ) t (cid:48) ) and p m,n to construct p m,n ( X ( m ) t , Y (2 n ) t , A ( m ) t (cid:48) + t − t (cid:48) , B ( n ) t (cid:48) + t − t (cid:48) , t ) = (cid:90) Λ +( m +2 n ) ˆ p ( X ( m ) t , Y (2 n ) t , t | X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) ) × S ,m ( t ; t (cid:48) , A ( m ) t (cid:48) ) S ,n ( t ; t (cid:48) , B ( m ) t (cid:48) ) p m,n ( X ( m ) t (cid:48) , Y (2 n ) t (cid:48) , A ( m ) t (cid:48) , B ( n ) t (cid:48) , t (cid:48) )d X ( m ) t (cid:48) d Y (2 n ) t (cid:48) . (9)Assuming that p m,n is continuous and differentiable, and the integration is interchangeable with differentiation inEq. (9), we take derivatives with respect to all variables t, X i , Y j , A i , B j to obtain ∂p m,n ∂t + m (cid:88) i =1 ∂ ( g ( X it , A it , t ) p m,n ) ∂X it + n (cid:88) j =1 ∂ ( g ( Y jt , B jt , t ) p m,n ) ∂Y jt + m (cid:88) i =1 ∂p m,n ∂A it + n (cid:88) j =1 ∂p m,n ∂B jt = − (cid:18) m (cid:88) i =1 β m,n ( A it , t ) + 2 n (cid:88) j =1 β m,n ( B jt , t ) (cid:19) p m,n + 12 m (cid:88) i =1 ∂ ( σ ( X it , A it , t ) p m,n )( ∂X it ) + 12 n (cid:88) j =1 ∂ ( σ ( Y jt , B jt , t ) p m,n )( ∂Y jt ) , (10)where p m,n ≡ p m,n ( X ( m ) t , Y (2 n ) t , A ( m ) t , B ( n ) t , t ). Hereafter, we will omit the subscript t for notational simplicity. Tofacilitate further analysis, we define a symmetrized density ρ m,n that is symmetric to the interchange of variables: ρ m,n ( X m , Y n , A m , B n , t ) = 12 n m ! n ! (cid:88) π n p m,n ( X ( m ) , π n ( Y (2 n ) ) , A ( m ) , B ( n ) , t ) (11)where A ( m ) = ( A ξ a (1) , . . . , A ξ a ( m ) ), B ( n ) = ( B ξ b (1) , . . . , B ξ b ( m ) ) are ordered ages, X ( m ) = ( X ξ a (1) , . . . , X ξ a ( m ) ), Y (2 n ) = ( Y ξ b (1) − , . . . , Y ξ b ( n ) ) are the corresponding sizes, and π n is some permutation Λ n → Λ n such that π n ( Y i ) , π n ( Y i − ) ∈ { Y i − , Y i } , π n ( Y i ) (cid:54) = π n ( Y i − ) , i = 1 , ..., n , i.e. , π n can interchange the sizes of twocells in a doublet. Therefore, there are 2 n total permutations π n . ξ a (1) , ..., ξ a ( m ) is a rearrangement such that A ξ a (1) ≥ A ξ a (2) ≥ ... ≥ A ξ a ( m ) and ξ b (1) , ..., ξ b ( n ) is a rearrangement such that B ξ b (1) ≥ B ξ b (2) ≥ ... ≥ B ξ b ( m ) .Defining such a ρ m,n allows us to remove the restriction that the ages must be presented in a descending order.Moreover, changing the order of two cells within in a doublet will not affect the value of ρ m,n . Definite integrals over ρ m,n are then related to those over p m,n via (cid:90) d X m d Y n d A m d B n ρ m,n ( X m , Y n , A m , B n , t ) = (cid:90) d X ( m ) d Y (2 n ) (cid:90) Λ d A ξ a (1) ...... (cid:90) A ξa ( m − d A ξ a ( m ) (cid:90) Λ d B ξ b (1) ... (cid:90) B ξb ( n − d B ξ b ( n ) p m,n ( X ( m ) , Y (2 n ) , A ( m ) , B ( n ) , t ) , (12)so ρ m,n is also a probability density distribution if p m,n is. Furthermore, the differential equation satisfied by ρ m,n for A m , B n > p m,n ∂ρ m,n ∂t + m (cid:88) i =1 ∂ρ m,n ∂A i + n (cid:88) j =1 ∂ρ m,n ∂B j + m (cid:88) i =1 ∂ ( ρ m,n g ( X i , A i , t )) ∂X i + n (cid:88) j =1 ∂ ( ρ m,n g ( Y j , B [ j +12 ] , t )) ∂Y j = − (cid:18) m (cid:88) i =1 β m,n ( A i , t ) + 2 n (cid:88) j =1 β m,n ( B j , t ) (cid:19) ρ m,n + 12 m (cid:88) i =1 ∂ ( σ ( X i , A i , t ) ρ m,n )( ∂X i ) + 12 n (cid:88) j =1 ∂ ( σ ( Y j , B [ j +12 ] , t ) ρ m,n )( ∂Y j ) . (13) B. Boundary Conditions
We now specify appropriate boundary conditions for ρ m,n that represent the birth of new cells with age zero. Byusing ordered ages, it is easy to derive the corresponding boundary conditions for p m,n defined in Eq. (9), which weomitted here, but which are nonzero if B n = 0 and zero if any entry in X m , Y n , A m , B k In this section, we will assume that ˜ β and β are independent of the population sizes m, n . Under this assumption,we are able to derive lower-dimensional ( e.g. , marginalized) projections of our kinetic theory (Eq. (13)) by averagingover a variable number of cell sizes: ρ ( h,k,(cid:96) ) m,n ( X h , A h , Y k +2 (cid:96) e , B k + (cid:96) , t ) = (cid:90) Λ d X h +1: m d Y o2 k +2 (cid:96) +1:2 n d A h +1: m d B k + (cid:96) +1: n ρ m,n , (16)where ρ m,n ≡ ρ m,n ( X m , Y n , A m , B n , t ), Λ ≡ Λ ( m − h )+(2 n − k − (cid:96) )+( m − h )+( n − k ) , and we define the notation X h +1: m :=( X h +1 , ..., X m ) , Y k +2 (cid:96) +1:2 n o := ( Y , Y , ..., Y k − , Y k +2 (cid:96) +1 , ..., Y n ), A h +1: m := ( A h +1 , ..., A m ), B k + (cid:96) +1: n :=( B k + (cid:96) +1 ,..., B n ) and Y k +2 (cid:96) e := ( Y , Y , ..., Y k , Y k +1 , Y k +2 , ..., Y k +2 (cid:96) ). The marginalized densities require three indices todescribe because although the size X m and age A m have a one-to-one correspondence for singlets, the twins, whilecarrying the same age, almost surely have different sizes due to asymmetric division and independent growth fluctu-ations immediately after birth. Thus, the number of ways to exit and enter each state depends on which types ofcells are “integrated over”. By marginalizing over Eq. (13), we find the kinetic equation satisfied by ρ ( h,k,(cid:96) ) m,n (in theremaining space X h , A h , Y k +2 (cid:96) e , B k > 0) becomes ∂ρ ( h,k,(cid:96) ) m,n ( X h , A h , Y k +2 l e , B k + (cid:96) , t ) ∂t + h (cid:88) i =1 ∂ρ ( h,k,l ) m,n ∂A i + k + (cid:96) (cid:88) j =1 ∂ρ ( h,k,(cid:96) ) m,n ∂B j + h (cid:88) i =1 ∂ ( g ( X i , A i , t ) ρ ( h,k,(cid:96) ) m,n ) ∂X i + k (cid:88) j =1 ∂ ( g ( Y j , A j , t ) ρ ( h,k,(cid:96) ) m,n ) ∂Y j + (cid:96) (cid:88) j =1 ∂ ( g ( Y k + j , A j , t ) ρ ( h,k,(cid:96) ) m,n ) ∂Y k + j − h (cid:88) i =1 ∂ ( σ ( X i , A i , t ) ρ ( h,k,(cid:96) ) m,n )( ∂X i ) − k (cid:88) j =1 ∂ ( σ ( Y j , B [ j +12 ] , t ) ρ ( h,k,(cid:96) ) m,n )( ∂Y j ) − (cid:96) (cid:88) j =1 ∂ ( σ ( Y k + j , B k +[ j +12 ] , t ) ρ ( h,k,(cid:96) ) m,n )( ∂Y k + j ) = − h (cid:88) i =1 β ( A i , t ) ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e , A h , B k + (cid:96) , t ) − k + (cid:96) (cid:88) j =1 β ( B j , t ) ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e , A h , B k + (cid:96) , t ) − ( m − h ) (cid:90) Λ d X h +1 d A h +1 β ( A h +1 , t ) ρ ( h +1 ,k,(cid:96) ) m,n ( X h +1 , Y k +2 (cid:96) e , A h +1 , B k + (cid:96) , t ) − n − k − (cid:96) ) (cid:90) Λ d Y k +2 d B k +1 β ( B k +1 , t ) ρ ( h,k +1 ,(cid:96) ) m,n ( X h , Y k +2 (cid:96) +2e , A h , B k + (cid:96) +1 , t )+ ( n − k − (cid:96) )( m + 1) n (cid:90) Λ d X h +1 d A h +1 β ( A h +1 , t ) ρ ( h +1 ,k,(cid:96) ) m +1 ,n − ( X h +1 , Y k +2 (cid:96) e , A h +1 , B k + (cid:96) , t )+ 2( n − k − (cid:96) )( m − h ) m (cid:90) Λ d Y k +2 d B k +1 β ( B k +1 , t ) ρ ( h,k +1 ,(cid:96) ) m − ,n ( X h , Y k +2 (cid:96) +2e , A h , B k + (cid:96) +1 , t )+ 2( n − k − (cid:96) ) m h (cid:88) i =1 β ( A i , t ) ρ ( h − ,k +1 ,(cid:96) ) m − ,n ( X h − i , Y k +2+2 (cid:96) e [ Y k +2 = X i ] , A h − i , B k + (cid:96) +1 [ B k +1 = A i ] , t ) , (17)and the associated boundary conditions become ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e [ Y k = y ] , A h , B k + (cid:96) [ B k = 0] , t ) = m + 1 n (cid:90) Λ d A h +1 d s ˜ β ( s + y, y, A h +1 , t ) ρ ( h +1 ,k − ,(cid:96) ) m +1 ,n − ( X h +1 [ X h +1 = s + y ] , Y k +2 (cid:96) − , A h +1 , B k + (cid:96) − , t )+ 2( m − h ) m (cid:90) Λ d B k d s ˜ β ( s + y, y, B k , t ) ρ ( h,k,(cid:96) ) m − ,n ( X h , Y k e [ Y k = s + y ] , A h , B k + (cid:96) , t )+ 2 m h (cid:88) i =1 (cid:90) Λ d s ˜ β ( s + y, y, A i , t ) ρ ( h − ,k − ,(cid:96) +1) m − ,n ( X h − i , Y k +2 (cid:96) e [ Y k +2 (cid:96) − = s + y, Y k +2 (cid:96) = X i ] , ...... A h − i , B k + (cid:96) [ B k = A i ] , t ) , (18) ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e [ Y k +2 (cid:96) − = y , Y k +2 (cid:96) = y ] , A h , B k + (cid:96) [ B k + (cid:96) = 0] , t ) = m + 1 n (cid:90) Λ d A h +1 ˜ β ( y + y , y , A h +1 , t ) ρ ( h +1 ,k,(cid:96) − m +1 ,n − ( X h +1 [ X h +1 = y + y ] , Y k +2 (cid:96) − , A h +1 , B k + (cid:96) − , t )+ 2( m − h ) m (cid:90) Λ d B k +1 ˜ β ( y + y , y , B k +1 , t ) ρ ( h,k +1 ,(cid:96) − m − ,n ( X h , Y k +2 (cid:96) e [ Y k +2 = y + y ] , A h , B k + (cid:96) , t )+ 2 m h (cid:88) i =1 ˜ β ( y + y , y , A i , t ) ρ ( h − ,k,(cid:96) ) m − ,n ( X h − i , Y k +2 (cid:96) e [ Y k +2 (cid:96) − = y + y , Y k +2 (cid:96) = X i ] , ...... A h − i , B k + (cid:96) [ B k + (cid:96) = A i ] , t ) , (19)and ρ ( h,k,(cid:96) ) m,n ( X h [ X i = 0] , Y k +2 (cid:96) e , A h , B k + (cid:96) , t ) = ρ ( h,k,(cid:96) ) m,n ( X h [ X i = ∞ ] , Y k +2 (cid:96) e , A h , B k + (cid:96) , t ) = 0 , i = 1 , , ..., h, (20) ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e [ Y j = 0] , A h , B k + (cid:96) , t ) = ρ ( h,k ) m,n ( X h , Y k +2 (cid:96) e [ Y j = ∞ ] , A h , B k + (cid:96) , t ) = 0 ,j = 2 , , ..., k, k + 1 , ..., k + 2 (cid:96), (21) ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e , A h [ A i = 0] , B k , t ) = 0 , i = 1 , , ..., h, (22) ρ ( h,k,(cid:96) ) m,n ( X h , Y k +2 (cid:96) e , A h , B k + (cid:96) , t ) = 0 , if two or more entries in B k + (cid:96) are 0 . (23)The first two terms on the RHS of Eq. (17) represent the division of a singlet/doublet in the current system whose ageis specified; the third and fourth terms on the RHS stand describe the division of a singlet and one cell of a doublet,respectively, whose age is not specified; the fifth term results from the division of a singlet, whose age and volume areunspecified, that induces the state transition ( m + 1 , n − → ( m, n ). The sixth term arises from division of one cellof a doublet that coverts the system from ( m − , n ) to ( m, n ). Finally, the last term represents the division of onecell in a doublet whose age is A i , ≤ i ≤ h and its undividing twin has size X i . In Eq. (18) and (19), the first termon their RHSs represent the division of a singlet, and the second term on their RHSs describe the division of one cellin a doublet, giving rise to a newborn doublet and leaving a singlet whose volume and age are integrated over. Thelast terms in the boundary conditions in Eq. (18) and (19) result from the division of a cell in a doublet, resulting ina newborn doublet and leaving a singlet whose volume and age are X i ∈ X h and A i ∈ A h , respectively.Our kinetic equations subsume all hierarchical equations for ρ ( h,k,(cid:96) ) m,n . First, we consider the lowest order equations( h = k = (cid:96) = 0) and the physical quantities that can easily be constructed such as the total number N = m + 2 n .The total expected cell population can be expressed as E [ N ( t )] = ∞ (cid:88) m =0 ∞ (cid:88) n =0 ( m + 2 n ) ρ (0 , , m,n , (24)which satisfiesd E [ N ( t )]d t = ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:20) m (cid:90) Λ d X d A β ( A , t ) ρ (1 , , m,n ( X , A , t ) + 2 n (cid:90) Λ d Y d B β ( B , t ) ρ (0 , , m,n ( Y , B , t ) (cid:21) = (cid:90) Λ d x d a β ( a, t ) n (1 , ( x, a, t ) (25)and involves the higher-dimensional densities ρ (1 , , m,n and ρ (0 , , m,n . The differential equation for E [ N ( t )] does not involvea any boundary condition, but it is not closed because it depends on n (1 , .Higher dimensional total number-density functions n ( k,(cid:96) ) ( x k , y (cid:96) , a k , b (cid:96) , t ) can also be generally defined: n ( k,(cid:96) ) = ∞ (cid:88) m =0 ∞ (cid:88) n =0 k (cid:88) r =0 (cid:88) ξ (0 ,r ) ∈ S k k + (cid:96) − r ( m ) r ( n ) k + (cid:96) − r ρ ( r,k − r,(cid:96) ) m,n ( X r [ X i = x ξ (0 ,r ) ( i ) ] , ......, Y k − r )+2 (cid:96) e [ Y j = x ξ ( r,k − r ) ( j ) , Y k − r )+ p = y p ] , A r [ A i = a ξ (0 ,r ) ( i ) ] , ......, B k − r + (cid:96) [ B j = a ξ ( r,k − r ) ( j ) , B k − r +[ p +12 ] = b [ p +12 ] ] , t ) , ≤ i ≤ r, ≤ j ≤ k − r, ≤ p ≤ (cid:96) (26)where x k := ( x , ..., x k ) , y (cid:96) := ( y , ..., y (cid:96) ) , a k := ( a , ..., a k ) , b (cid:96) := ( b , ..., b (cid:96) ), ( m ) r = m ! / ( m − r )! is the fallingfactorial, S k = { , , ..., k } . The sum (cid:80) ξ (0 ,r ) ∈ S k includes summing over all elements ξ (0 ,r ) ∈ Ω r , the set that containsall possible choices of choosing r elements in S k , and ξ ( r,k − r ) := ( ξ ( r + 1) , ξ ( r + 2) , ...ξ ( k )) = S k \ ξ (0 ,r ) . We require ξ (0 ,r ) ( i ) < ξ (0 ,r ) ( j ) , ξ ( r,k − r ) ( i ) < ξ ( r,k − r ) ( j ) , ∀ i < j and r ≤ m, k − r ≤ n in Eq. (26). With a β independent of m, n ,the PDE satisfied by n ( k,(cid:96) ) ( x k , y (cid:96) , a k , b (cid:96) , t ) is ∂n ( k,(cid:96) ) ∂t + k (cid:88) i =1 ∂n ( k,(cid:96) ) ∂a i + (cid:96) (cid:88) j =1 ∂n ( k,(cid:96) ) ∂b j + k (cid:88) i =1 ∂ ( n ( k,(cid:96) ) g ( x i , a i , t )) ∂x i + (cid:96) (cid:88) j =1 ∂ ( n ( k,(cid:96) ) g ( y j , b [ j +12 ] , t )) ∂y j = − (cid:18) k (cid:88) i =1 β ( a i , t ) n ( k,(cid:96) ) + (cid:96) (cid:88) j =1 β ( b j , t ) (cid:19) n ( k,(cid:96) ) + 12 k (cid:88) i =1 ∂ ( n ( k,(cid:96) ) σ ( x i , a i , t ))( ∂x i ) + 12 (cid:96) (cid:88) j =1 ∂ ( n ( k,(cid:96) ) σ ( y j , b [ j +12 ] , t ))( ∂y j ) , (27)along with the boundary conditions n ( k,(cid:96) ) ( x k [ x v = x ] , a k [ a v = 0] , y (cid:96) , b (cid:96) , t ) = ∞ (cid:88) m =0 ∞ (cid:88) n =0 k − (cid:88) r =0 (cid:88) ξ (0 ,r ) ∈ S − vk (cid:96) + k − r ( m ) r ( n ) k + (cid:96) − r × ρ ( r,k − r,(cid:96) ) m,n ( X r [ X i = x ξ (0 ,r ) ( i ) ] , Y k − r +2 (cid:96) e [ Y j = x ξ ( r,k − r ) ( j ) , Y k + p = y p ] , ......, A r [ A i = a ξ (0 ,r ) ( i ) ] , B (cid:96) + k − r [ B j = a ξ ( r,k − r ) ( j ) , B k − r +[ p +12 ] = b [ p +12 ] ] , t )= 2 (cid:90) Λ d s d a ˜ β ( x + s, x, a, t ) n ( k,(cid:96) ) ( x k [ x k = x + s ] , y (cid:96) , a k [ a k = a ] , b (cid:96) , t )+ 2 k (cid:88) u =1 , (cid:54) = v (cid:90) Λ d s ˜ β ( x + s, x, a u , t ) n ( k − ,(cid:96) +1) ( x k − u, − v , a k − u, − v , ......, y (cid:96) +2 [ y (cid:96) +1 = x u , y (cid:96) +2 = s + x ] , b (cid:96) +1 [ b (cid:96) +1 = a u ] , t ) (28) n ( k,(cid:96) ) ( x k , y (cid:96) [ y v − = y , y v = y ] , a k , b (cid:96) [ b v = 0] , t ) = ∞ (cid:88) m =0 ∞ (cid:88) n =0 k (cid:88) r =0 (cid:88) ξ (0 ,r ) ∈ S k (cid:96) + k − r ( m ) r ( n ) k + (cid:96) − r × ρ ( r,k − r,(cid:96) ) m,n ( X r [ X i = x ξ (0 ,r ) ( i ) ] , Y (cid:96) + k − r e [ Y j = x ξ ( r,k − r ) ( j ) , Y k + q = y q ] , ......, A r [ A i = a ξ (0 ,r ) ( i ) ] , B (cid:96) + k − r [ B j = a ξ ( r,k − r ) ( j ) , B k − r +[ q +12 ] = b [ q +12 ] ] , t )= 2 (cid:90) Λ d a ˜ β ( y + y , y , a, t ) n ( k +1 ,(cid:96) − ( x k +1 [ x k +1 = y + y ] , y (cid:96) − (2 v − , − v , a (cid:96) +1 [ a (cid:96) +1 = a ] , b (cid:96) − v , t )+ 2 k (cid:88) u =1 , (cid:54) = v ˜ β ( y + y , y , a u , t ) n ( k − ,(cid:96) ) ( x k − u , a k − u , y (cid:96) [ y v − = y + y , y v = x u ] , b (cid:96) [ b v = a u ] , t ) , (29)where x k − u := ( x , ..., x u − , x u +1 , ..., x k ), a k − u := ( a , ..., a u − , ..., a u +1 , ...a k ), x k − u, − v := ( x , ..., x u − , x u +1 , ..., x v − , x v +1 ,..., x k ), a k − u, − v := ( a , ..., a u − , a u +1 , ..., a v − , a v +1 , ..., a k ), y (cid:96) − (2 v − , − v := ( y , ..., y v − , y v +1 , ..., y (cid:96) ), b (cid:96) − v :=( b , ..., b v − ,b v +1 , ..., b (cid:96) ) and S − vk := { , , ..., v − , v + 1 , .., k } . The additional conditions, n ( k,(cid:96) ) ( x k , a k , y (cid:96) , b (cid:96) , t ) = 0 (cid:26) if any x i , y j = 0 , ∞ if two or more a i or b j =0 (30)are found by using Eq. (26) in Eqs. (19). Note that if we take k = 1 , (cid:96) = 0, with an m, n -independent β , the “1-point”total mean population density n (1 , ( x, a, t ) in volume x and age a at time t is simply n (1 , ( x, a, t ) ≡ ∞ (cid:88) m =1 ∞ (cid:88) n =0 mρ (1 , , m,n ( X [ X = x ] , A [ A = a ] , t ) + ∞ (cid:88) m =0 ∞ (cid:88) n =1 nρ (0 , , m,n ( Y [ Y = x ] , B [ B = a ] , t ) , (31)and obeys a first-moment (in both dimension and particle number), closed PDE ∂n (1 , ∂t + ∂n (1 , ∂a + ∂ ( gn (1 , ) ∂x = − β ( a, t ) n (1 , ( x, a, t ) + 12 ∂ ( σ n (1 , ) ∂x (32)with associated boundary conditions specified at a = 0 , x = 0 , x = ∞ n (1 , ( x, , t ) = 2 n ∞ (cid:88) m =0 ∞ (cid:88) n =1 ρ (0 , , m,n ( Y [ Y = x ] , B [ B = 0] , t )= 2( m + 1) ∞ (cid:88) m =0 ∞ (cid:88) n =1 (cid:90) ∞ x d X (cid:90) Λ d A ˜ β ( X , x, A , t ) ρ (1 , , m +1 ,n − ( X , A , t )+ 4 n ∞ (cid:88) m =0 ∞ (cid:88) n =1 (cid:90) ∞ x d Y (cid:90) Λ d B ˜ β ( Y , x, B , t ) ρ (0 , , m,n ( Y , B , t )= 2 (cid:90) ∞ x d z (cid:90) Λ d a ˜ β ( z, x, a, t ) n (1 , ( z, a, t ) ,n (1 , (0 , a, t ) = n (1 , ( ∞ , a, t ) = 0 . (33) (a) (b) (c) FIG. 1: A map of boundary condition interdependences for single-density kinetic theory. In (a) we indicate the dependenceof the boundary condition for n ( k,(cid:96) ) ( x k , a k , y (cid:96) , b (cid:96) , t ) if any a i = 0. The boundary condition for n ( k,(cid:96) ) depends on itself and n ( k − ,(cid:96) +1) ; for example, n (0 , is required for the boundary condition for n (2 , , so the red arrow points from n (0 , to n (2 , .In (b) we indicate the dependence of the boundary condition for n ( k,(cid:96) ) ( x k , a k , y (cid:96) , b (cid:96) , t ) if any b j = 0. Here, the boundarycondition for n ( k,(cid:96) ) depends on n ( k +1 ,(cid:96) − and n ( k − ,(cid:96) ) . (c) An example of an explicit sequence of calculations to find n (1 , starting from n (1 , . Note that the PDEs for all multi-point single-density functions n ( k,(cid:96) ) are closed. However, the boundary conditionscouple n ( k,(cid:96) ) , k + (cid:96) > n ( k +1 ,(cid:96) − , n ( k − ,(cid:96) ) , or n ( k − ,(cid:96) +1) . Thus, although the full models for n ( k,(cid:96) ) , k + (cid:96) > n ( k (cid:48) ,(cid:96) (cid:48) ) such that k (cid:48) +2 (cid:96) (cid:48) ≤ k +2 (cid:96) , and therefore all n ( k,(cid:96) ) , k + (cid:96) > n (1 , . For instance, we can calculate n (0 , from n (1 , , and then n (2 , , n (1 , , and so on. How the different n ( k,(cid:96) ) are connected through the boundary conditions are illustrated in Fig. 1,demonstrating the sequence to follow to fully solve the single-density problem. The differential equation satisfied bythe lowest order moment E [ N ( t )] requires n (1 , , as indicated by the shaded blue arrow in Fig. 1(a). In Fig. 1(c) weshow a sequence of boundary condition calculations to find n (1 , : the equations satisfied by n (1 , are fully closedso n (1 , can be first calculated. In the second step, we use n (1 , to construct the boundary condition and solve for n (0 , . The third step is to use n (0 , to construct the boundary condition and solve for n (2 , . The boundary conditiondependence of n (1 , , n (2 , is indicated by blue arrows. The forth step and fifth steps are to solve for n (1 , and n (3 , ,whose boundary condition dependences are indicated by the green arrows. Next, we calculate n (2 , , n (0 , , and finally n (1 , , whose boundary condition dependences are shown by the red arrows.These higher dimensional results capture the stochasticity arising only from noisy growth of each cell (through thediffusive terms in Eqs. (27) and (32)). The demographic stochasticity arising from random birth (and death) times0affects the total population and is most directly probed by higher number correlations. For example, the differentialequation satisfied by E [ N ( t )] = ∞ (cid:88) m =0 ∞ (cid:88) n =0 ( m + 2 n ) ρ (0 , , m,n , (34)is d E [ N ( t )]d t = ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:20) (2 m + 4 mn + m ) (cid:90) d X d A β ( A , t ) ρ (1 , , m,n ( X , A , t )+ (8 n + 4 mn + 2 n ) (cid:90) d Y d B β ( B , t ) ρ (0 , , m,n ( Y , B , t ) (cid:21) . (35)The lowest order Eq. (25) decouples for β ( t ) which does not depend on age, trivially reducing to d E [ N ( t )] / d t = β ( t ) E [ N ( t )]. As for E [ N ( t )], if the division rate function is not dependent on age, Eq. (35) reduces tod E [ N ( t )]d t = 2 β ( t ) E [ N ( t )] + β ( t ) E [ N ( t )] . (36)It is also possible to derive the differential equations satisfied by any d E [ N k ( t )] / d t, k ∈ N + starting from Eq. (17).Such equations, as well as those for higher number-moments such as (cid:80) m,n m k ρ ( h,k,(cid:96) ) m,n are not closed and form complexhierarchies that need additional assumptions to close. IV. GENERALIZATIONSA. Incorporation of death Here, we show how our kinetic theory is modified when an age and size-dependent death, occuring with rate µ ( a, t ),is incorporated. By defining γ ( a, t ) = β ( a, t ) + µ ( a, t ) (37)the joint survival probabilities S ,m and S ,n in Eq. (6) are modified by˜ S ,m ( t ; t (cid:48) , A mt (cid:48) ) = m (cid:89) i =1 e − (cid:82) tt (cid:48) γ ( A it (cid:48) − t (cid:48) + s,s )ds , ˜ S ,n ( t ; t (cid:48) , B nt (cid:48) ) = n (cid:89) j =1 (cid:104) e − (cid:82) tt (cid:48) γ ( B jt (cid:48) − t (cid:48) + s,s )ds (cid:105) . (38)Following the previous derivations, we find ∂ρ m,n ∂t + m (cid:88) i =1 ∂ρ m,n ∂A i + n (cid:88) j =1 ∂ρ m,n ∂B j + m (cid:88) i =1 ∂ ( g ( X i , A i , t ) ρ m,n ) ∂X i + n (cid:88) j =1 ∂ ( g ( Y j , B j , t ) ρ m,n ) ∂Y j = − (cid:18) m (cid:88) i =1 γ ( A i , t ) + 2 n (cid:88) j =1 γ ( B j , t ) (cid:19) ρ m,n + 12 m (cid:88) i =1 ∂ ( σ ( X i , A i , t ) ρ m,n )( ∂X i ) + 12 n (cid:88) j =1 ∂ ( σ ( Y j , B j , t ) ρ m,n )( ∂Y j ) + ( m + 1) (cid:90) Λ d A m +1 d X m +1 µ ( A m +1 , t ) ρ m +1 ,n ( X m +1 , Y n , A m +1 , B n , t ) (39)+ 2( n + 1) m m (cid:88) i =1 (cid:90) Λ d x µ ( A i , t ) ρ m − ,n +1 ( X m − i , Y n +2 [ Y n +1 = x, Y n +2 = X i ] , A m − i , B n +1 [ B n +1 = A i ] , t ) , where the argument of ρ m,n in the first two lines is ( X m , Y n , A m , B n , t ).1The boundary conditions for ρ m,n are the same as Eq. (14) and Eq. (15) since only cell division contributes to theboundary term, and no cell can have 0 or infinitely large volume at any time. Similarly, we can define the marginaldistribution ρ ( h,k,l ) m,n ( X h , Y k +2 l e , A h , B k , t ) and the population density function with respect to volume x and age a attime t is n (1 , ( x, a, t ) = ∞ (cid:88) m =1 ∞ (cid:88) n =0 mρ (1 , , m,n ( X [ X = x ] , A [ A = a ] , t ) + ∞ (cid:88) m =0 ∞ (cid:88) n =1 nρ (0 , , m,n ( Y [ Y = x ] , B [ B = a ] , t ) . (40)By similar calculations as in Section (III), we obtain the differential equation satisfied by n (1 , ( x, a, t ) ∂n (1 , ∂t + ∂ ( gn (1 , ) ∂x + ∂n (1 , ∂a − ∂ ( σ n (1 , )( ∂x ) = − ( β ( a, t ) + µ ( a, t )) n (1 , ( x, a, t ) , (41)with boundary conditions specified at a = 0 and x = 0 , ∞ n (1 , ( x, , t ) = 2 n ∞ (cid:88) m =0 ∞ (cid:88) n =1 ρ (0 , , m,n ( Y [ Y = x ] , B [ B = 0] , t )= 2( m + 1) ∞ (cid:88) m =0 ∞ (cid:88) n =1 (cid:90) ∞ x d X (cid:90) Λ d A ˜ β ( X , x, A , t ) ρ (1 , , m +1 ,n − ( X , A , t )+ 4 n ∞ (cid:88) m =0 ∞ (cid:88) n =1 (cid:90) ∞ x d Y (cid:90) Λ d B ˜ β ( Y , x, B , t ) ρ (0 , , m − ,n ( Y , B , t )= 2 (cid:90) ∞ d a (cid:90) ∞ x d z ˜ β ( z, x, a, t ) n (1 , ( z, a, t ) ,n (1 , (0 , a, t ) = n (1 , ( ∞ , a, t ) = 0 . (42) B. Correlated noise in growth rate In this subsection we consider a model in which the noise in growth rates are correlated across cells. By defining Z m, n = ( X m , Y n ) and C m, n = ( A m , B , B , ..., B n , B n ) to be the volumes and ages of m singlets and n doubletsat time t , we can describe the growth rate asd Z m, nt = G m, n ( Z m, nt , C m, nt , t )d t + Σ m, n ( Z m, nt , C m, nt , t )d W pt , (43)where G m, n ∈ R m +2 n , Σ m, n ( Z m, nt , C m, nt , t ) = ( σ ) ij ∈ R ( m +2 n ) × p and W pt is a p -dimensional i.i.d standard Wienerprocess [22]. For simplicity, we assume that the i th component of G m, n is g i ( Z it , C it , t ) = g ( Z i , C i , t ), indicatingthat the deterministic part of the growth rate is identical for all cells. We further assume that the variance ingrowth rates for all cells is identical: (cid:80) p(cid:96) =1 σ i,(cid:96) = σ , ∀ i . Following our derivation in Section (II), we find that ρ m,n ( X m , Y n , A m , B n , t ) satisfies ∂ρ m,n ∂t + m (cid:88) i =1 ∂ρ m,n ∂A i + n (cid:88) j =1 ∂ρ m,n ∂B i + m (cid:88) i =1 ∂ ( g ( t, X i , A i ) ρ m,n ) ∂X i + n (cid:88) j =1 ∂ ( g ( t, Y j , B [( j +1) / ) ρ m,n ) ∂Y j = − (cid:18) m (cid:88) i =1 β ( A i , t ) + n (cid:88) j =1 β ( B j , t ) (cid:19) ρ m,n ( X m , Y n , A m , B n , t ) + m +2 n (cid:88) s ,s =1 ∂ ( ρ m,n D s ,s ) ∂Z s ∂Z s , (44)where D s ,s = (cid:80) p(cid:96) =1 σ s ,(cid:96) σ s ,(cid:96) . The boundary conditions for ρ m,n are the same as that described by Eq. (14) andEq. (15). Similarly, we can define the marginal distribution density function ρ ( h,k,(cid:96) ) m,n in the same way as in Section 3,and it can be verified that the differential equations as well as the boundary conditions satisfied by ρ (1 , , m,n ( X [ X = x ] , A [ A = a ] , t ) , ρ (0 , , m,n ( X [ X = x ] , A [ A = a ] , t ) are the same as those satisfied by ρ (1 , , m,n ( X [ X = x ] , A [ A = a ] , t ) and ρ (0 , , m,n ( Y [ Y = x ] , B [ B = a ] , t ) in Eq. (17) and Eq. (19), although the differential equations satisfied2by ρ m,n in Eq. (44) and in Eq. (11) are different. The equation and boundary conditions for the “1-point” densityfunction n (1 , ( x, a, t ) are identical to those in Eq. (32) and Eqs. (33) since correlations are not captured by a mean-field description of only one coordinate ( x, a ). The differences between correlated and uncorrelated growth noiseamong cells may arise in the differential equations for n ( k,(cid:96) ) ( x k , a k , y (cid:96) , b (cid:96) , t ) , (cid:96) + k ≥ V. SUMMARY AND CONCLUSIONS In this paper, we rigorously constructed a kinetic theory for structured populations, in particular for age- and size-structured cell proliferation models. We considered stochasticity in both an individual cell’s growth rate (“intrinsic”stochasticity) and the cell number fluctuations from random birth and death event times (“demographic” stochastic-ity). Derivations of the kinetic theory requires separation of ’singlet’ and ’doublet’ populations, as was proposed in[16]. However, taking into account both the size and age dependence as well as randomness in growth rates leads tothe much more complex computation which we performed here.One of our main results are the kinetic equations and boundary conditions described by Eqs. (13), (14), and (15).Marginalized densities are also found to obey more complex equations that form a hierarchy (Eqs. (17), (19), and(23)). By taking single-density averages over these equations, we find closed PDEs that govern multi-point densityfunctions (Eq. (27)). However, the associated boundary conditions, Eq. (28), couple density functions of differentdimensions. Nonetheless, density function of all dimensions can be successively solved starting from the “1-point”density n (1 , ( x, a, t ) which obeys Eqs. (32) and (33), a 2+1-dimensional second order PDE and boundary conditionthat is analogous to the classic McKendrick equation but that a includes a diffusive size term arising from stochasticityin growth rates. The explicit equations for the first and second moments of the total population are given by Eqs. (25)and (35), respectively.Generalizations and extensions to our basic kinetic theory are also investigated. For example, we derived the kineticequations when a Markovian age-dependent death process is included (Eqs. (39), and (41), (42)). We also considerednoise in growth rates that are correlated across cells and showed these effects arising in “cross-diffusion” terms in theassociated kinetic (and higher moment) equations.Our unifying kinetic theory enables one to systematically analyze cell populations at both the individual andpopulation levels. A full kinetic theory may be useful for studying other processes such as failure in multicomponentsystems that age and evolve [23]. Further extensions of our kinetic equations that are feasible are to include spatialdistribution [24] or correlations in growth rates across generations [13]. It is also possible to consider stochasticityfor different cell division strategies [21]. Finally, efficient numerical methods for solving our kinetic equations can bedeveloped, for instance in [25] equations similar to Eq. (32) and Eq. (33) which describes the dynamics of n (1 , aresolved accurately and efficiently. Acknowledgements This research was made possible through funding support from the Army Research Office (W911NF-18-1-0345),the NIH (R01HL146552), and the National Science Foundation (DMS-1814364). References [1] von Foerster H 1959 The Kinetics of Cellular Proliferation, Grune and Stratton Current Biology Bulletin of the American Physical Society [4] Robert L, Hoffmann M, Krell N, Aymerich S, Robert J and Doumic M 2014 BMC Biology The Dynamics of Physiologically Structured Populations (Springer)[7] Sompayrac L and Maaloe O 1973 Nature: New Biology Nature Current Biology PLoS ONE e0182633 [11] Wessels J G H 1994 Annual Review of Phytopathology Biophysical Journal SIAM Journal on Applied Mathematics Kinetic and Related Models Physical Review E Journal of Statistical Physics Journal of Statistical Mechanics IEEE Life Sciences Letters Annual Review of Biophysics Physical Review E (4) 042139[21] Nieto C, Vargas-Garcia C and Pedraza J M 2020 bioRxiv:2020.09.29.319251 [22] Durrett R 2005 Cambridge U Press Journal of The Royal Society Interface Structured Population Models in Biology and Epidemiology vol 1936 (Springer)[25] Xia M, Shao S and Chou T 2020 arXiv:2009.13170 Appendix: conservation of probability We now define probability fluxes J m,n ; m +1 ,n − ( t ) = ( m + 1) (cid:90) d X m d Y n − d A m d B n − (cid:90) Λ d y d y d s ˜ β m +1 ,n − ( y + y , y , s, t ) × ρ m +1 ,n − ( X m +1 [ X m +1 = y + y ] , Y n − , A m +1 [ A m +1 = s ] , B n − , t ) ,J m,n ; m − ,n ( t ) = 2 nm (cid:90) d X m d Y n − d A m d B n − (cid:90) Λ d y d y m (cid:88) i =1 ˜ β m − ,n ( y + y , y , A i , t ) × ρ m − ,n ( t, X m − i , Y n [ Y n − = X i , Y n = y + y ] , A m − i , B n [ B n = A i ] , t ) ,J m,n ; m (cid:48) ,n (cid:48) ( t ) = 0 , if m + 2 n − m (cid:48) − n (cid:48) (cid:54) = 1 . (45) J m,n ; m (cid:48) ,n (cid:48) ( t )d t is the probability flux within time [ t, t + d t ] from state ( m (cid:48) , n (cid:48) ) to state ( m, n ) arising from from celldivision. When d t is sufficiently small, the probability that more than one cell divides during [ t, t + d t ] is o (d t ),which is negligible, allowing us to set J m,n ; m (cid:48) ,n (cid:48) ( t ) = 0 if m + 2 n − m (cid:48) − n (cid:48) (cid:54) = 1. We now verify the conservation ofprobability flux J m − ,n +1; m,n ( t ) + J m +1 ,n ; m,n ( t )= (cid:90) d X m d Y n d A m d B n (cid:18) m (cid:88) i =1 β m,n ( A i , t ) ρ m,n + n (cid:88) i = j β m,n ( B j , t ) ρ m,n = (cid:90) d X m d Y n d A m d B n (cid:18) mβ m,n ( A m , t ) ρ m,n + 2 nβ m,n ( B n , t ) ρ m,n (cid:19) , (46)where ρ m,n = ρ m,n ( X m , Y n , A m , B n , t ). The first term is J m − ,n +1; m,n ( t ) = m (cid:90) d X m − d Y m d A m − d B n (cid:90) Λ d y d y d A m ˜ β m,n ( y + y , y , A m , t ) × ρ m,n ( X m [ X m = y + y ] , Y n , A m , B n , t )= m (cid:90) d X m − d Y n d A m − d B n (cid:90) Λ d A m d( y + y ) (cid:90) y + y d y ˜ β m,n ( y + y , y , A m , t ) × ρ m,n ( X m [ X m = y + y ] , Y n , A m , B n , t )= m (cid:90) d X m − d Y n d A m − d B n (cid:90) Λ d A m d X m β m,n ( A m , t ) ρ m,n (47)4which is exactly the first term on the right hand side of Eq. (46). The second term J m +1 ,n ; m,n ( t ) = 2 nm + 1 (cid:90) d X m +1 d Y n d A m +1 d B n − (cid:90) Λ d y d y m +1 (cid:88) i =1 ˜ β m,n ( y + y , y , A i , t ) × ρ m,n ( X m +1 − i , Y n [ Y n − = X i , Y n = y + y ] , A m +1 − i , B n [ B n = A i ] , t )= 2 nm + 1 m +1 (cid:88) i =1 (cid:90) d X m +1 d Y n − d A m +1 d B n − (cid:90) Λ d( y + y ) (cid:90) y + y d y ˜ β m,n ( y + y , y , A i , t ) × ρ m,n ( X m +1 − i , Y n [ Y n − = X i , Y n = y + y ] , A m +1 − i , B n [ B n = A i ] , t )= 2 n (cid:90) d X m d Y n d A m d B n β m,n ( B n , t ) ρ m,n (48)which is precisely the second term on the right hand side of Eq. (46). We have thus verified that the probabilityflux out of state ( m, n ) due to cell division is the sum of probability currents into ( m − , n + 1) and into ( m + 1 , n ).Summing up over m and n , we obtain for m + n > ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:18) J m − ,n +1; m,n ( t ) + J m +1 ,n ; m,n ( t ) (cid:19) = ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:90) d X m d Y n d A m d B n (cid:18) mβ m,n ( A m , t ) ρ m,n + 2 nβ m,n ( B n , t ) ρ m,n (cid:19) . (49)Finally, it is readily observed that ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:90) d X m d Y n d A m d B n ∂ρ m,n ∂t = ∞ (cid:88) m =0 ∞ (cid:88) n =0 n (cid:88) j =1 (cid:90) d X m d Y n d A m d B n − j ρ m,n ( X m , Y n , A m , B n [ B j = 0] , t ) − ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:90) d X m d Y n d A m d B n (cid:18) m (cid:88) i =1 β m,n ( A i , t ) ρ m,n + n (cid:88) j =1 β m,n ( B j , t ) ρ m,n (cid:19) = ∞ (cid:88) m =1 ∞ (cid:88) n =0 ( J m,n ; m − ,n − J m − ,n +1; m,n ) − ∞ (cid:88) m =0 ∞ (cid:88) n =0 J m +1 ,n ; m,n + ∞ (cid:88) m =0 ∞ (cid:88) n =1 J m,n ; m +1 ,n − = 0 (50)Therefore, we have verified that ∞ (cid:88) m =0 ∞ (cid:88) n =0 (cid:90) d X m d Y n d A m d B n ρ m,nm,n