Numerical methods to solve age structured equations A tutorial for biologists E. Grenier
Introduction: what is an age structured model ? In some cases it is interesting to detail the repartition of a population in age classes. For instance: • Structure in age of human population: so called age pyramid. • Patients may be classed according to the duration of their illness • Populations of cells which undergo apoptosis may be classed according to their degree of avancement in this process. Therefore • Populations may be classed according to their age • Populations undergoing a process may be classed according to its degree of completion. which happens often !
Age pyramid : introduction 1.1 Introduction Let us consider a population of individues, of various ages, for instance men. The possible ages vary from 0 to say 100 years. Let f(t,a) denotes the density of men which at time t, precisely have age a, that is, if a and b are close enough, the number of men of age between a and b is nearly f(t,a) (b-a). Note that here a is a continuous variable, and not a discrete one: instead of simply considering age as being 1, 2, 3, …. years, we consider age like 1,234343… year and so on. This turns to be more convenient when making computations. Let Δt be a small time interval. Men of age a at time t are of age a+Δt at time t + Δt, except if they die meanwhile. Let m(a) be the mortality corresponding to age a. Then f(t+dt,a+dt) = f(t,a) – f(t,a) m(a) Δt which leads to ∂t f(t,a) + ∂a f(t,a) = - m(a) f(t,a)
Age pyramid : introduction This equation must be completed by boundary conditions. Namely at time t, f(t,0) is simply the number of newborn babies, which we will call b(t). It must also be completed by the initial data f(0,a). ∂t f(t,a) + ∂a f(t,a) = - m(a) f(t,a) f(t,0) = b(t)
Age pyramid : numerical algorithms 2. Numerical algorithms We discretize • Time: let Δt be the time step. We compute approximate values of f at times t_n = n Δt • Age: let Δa be the age step. We compute approximate values for ages a_j = j Δa We then discretize the equation in ( f(t_n + Δt,a_j) – f(t_n,a_j) )/ Δt + ( f(t_n, a_j) – f(t_n, a_j - Δa) )/ Δa = - m(a_j) f(t_n,a_j) with f(t_n,a_0) = b(t_n).
Age pyramid : numerical algorithms 3. Matlab program amax = 100; % maximal age da = 0.2; % age step N = amax/da; % number of age classes tmax = 200; % maximal time dt = 0.1; % time step mort = 0.1; % mortality F = zeros(1,N); % initialization of f for J=1:50/da, F(J) = 1; % initial data end for t=dt:dt:tmax, F(2:N) = F(2:N)*(1-dt/da) + F(1:N-1)*dt/da; % numerical algorithm F(1:N) = F(1:N)-mort*F(1:N).*([1:N] > N/2)*dt; % mortality F(1) = 1 + sin(t/5); % natality plot([0:N-1]*da + i*F); % plot drawnow; end
Age pyramid : a movie 4. Illustration Initial data : uniform repartition between 0 and 50 years Boundary data: sinusoidal: 1 + sin(t/5) Click to run the movie
Age pyramid : comments 5. Comments The numerical scheme is a litte diffusive
Age structured models • Introduction Let us consider a population of ill animals. We describe the degree of this illness by a, ranging from 0 (beginning of the illness) to 1 (end of the illness). The progression in the illness is not necessary uniform, and may depend on a. Let c(a) be the speed of progression of the illness. The case where c(a) = 1 corresponds to the age pyramid. The case where c depends on a is a little different. Let us denote by f(t,a) the density of animals in the state a.
Age structured models: equation 2. Equation Let us consider a small interval [a,a+da] and let us count the animals in this interval at time t and t+dt where da and dt are small. At time t: the number of animals in [a,a+da] is f(t,a) da At time t+dt, it is f(t+dt,a) da Meanwhile, some animals reach this state of illness. All the animals in [a-c(a)dt,a] enter this interval during time dt, therefore f(t,a) c(a) dt Animals reach the state a. Similarly f(t,a+da) c(a+da) dt get into a state worse than a+da. Therefore f(t+dt,a) da – f(t,a) da = f(t,a) c(a) dt – f(t,a+da) c(a+da) dt
Age structured models: equation This leads to ∂t f + ∂a (c f) = 0. This is discretized in ( f(t_n + dt, a_j) – f(t_n, a_j) ) / dt + ( c(a_j) f(t_n, a_j) – c(a_(j -1)) f(t_n, a_(j – 1)) / dx = 0 where f(t_n,0) is given (number of animals which get ill at time t_n).
Age structured models: Matlab program 3. Matlab program da = 0.02; % age step N = 1/da; % number of age classes tmax = 2; % maximal time of simulation dt = 0.01; % time step c = zeros(1,N); % initialisation of c c = 1 - [0:da:1-da].^2; % c= 1 – a^2 F = zeros(1,N); % initialisation of the initial data for J=1:N/2, % definition of initial data F(J) = 1; end for t=dt:dt:tmax, F(2:N) = F(2:N).*(1-c(2:N)*dt/da) + F(1:N-1).*c(1:N-1)*dt/da; % algorithm F(1) = 1; % data at a = 0 plot([0:N-1]*da + i*F); % plot drawnow; end
Age structured models: Matlab program Click to launch. Vertical scale changes with time: as c vanishes for a = 1, animals get more and more ill but never reach a = 1. Therefore they accumulate near a = 1.