first , we prove that both markov chains defined by ( 2 ) , ( 3 ) are ergodic and , second , we derive formula ( 4 ) . the ergodicity is a direct consequence of foster - lyapunov criteria ( corollary 8.7 p. 214 in robert ( 2003 ) or proposition 1.3 in hairer ( 2016 ) , for instance ) . the chains are time - homogeneous and , for l0>a : e(l1l0)=e(ab0+c0g0)=a/2+f(l0)e(g0),and this last quantity is negative for large enough values of l0 because f tends to zero , either f is f(l)=(1+(lls))1 or f(l)=1lis . this proof of ergodicity works for any function f having a limit l at infinity which satisfies l < moreover , the markov chains ( ln ) defined by ( 2 ) , ( 3 ) are irreducible and aperiodic . therefore , in both cases , there exists a unique equilibrium distribution , denoted l. for the model ( 3 ) , it is characterized by the fact that l and ( lab0)++g01{lis } have the same probability generating function . to establish formula ( 4 ) we distinguish the three regimes : l<ab0 , ab0lis , l>is , and the cases b0=0 or b0=1 . for u such that |u|1 , we obtain:(a.1)2e(ul)=2e(u(lab0)++g01lis)=2e(ug01{l<ab0})+2e(ulab0+g01{ab0lis})+2e(ulab01{l>is})=e(ug0)[p(l<a)+p(l<0)]+e(ug0)[e(ula1{alis})+e(ul1{0lis})]+[e(ula1{l>is})+e(ul1{l>is})]=e(ug0)sa1(1)+e(ug0)[ua(sis(u)sa1(u))+sis(u)]+[ua(e(ul)sis(u))+(e(ul)sis(u ) ) ] . from ( a.1 ) , we get : e(ul)(1ua)=sis(u)(e(ug0)1)(1+ua)+e(ug0)(sa1(1)uasa1(u ) ) . as the probability generating function of a geometric distribution is explicitly given by e(ug0)=p(1uq)1 , with q=1p , we obtain ( 4 ) after multiplication by ( 1uq)ua . in order to compute the k for all k , we identify the coefficients of the power series of each side of ( 4 ) distinguishing cases for the values of k. for simplicity we set k=0 for k<0 . using the identities : ( ua1)(1uq)=qua+1+ua+qu1 , 
 ( u1)(ua+1)=ua+1ua+u1 , we get:(a.2)k=0(qka1+ka+qk1k)uk=k=0k(quk+a+1+uk+a+quk+1uk)=(qua+1+ua+qu1)e(ul)=q(ua+1ua+u1)sis(u)+psa1(1)uapsa1(u)=q[k = a+1is+a+1ka1ukk = ais+akauk+k=1is+1k1ukk=0iskuk]+psa1(1)uapk=0a1kuk . the following table gives the recurrence relations obtained after identification in the coefficients of ( a.2).tablekrelation0,a1aa=21ppoa1k=1ka+1,isk=21ppka1 + 2ppkais+1is+1=2qisa+(1+q)is+1ais+2,is+ak=2qka1+(1+q)ka+qk1is+a+1is+a+1=2qis+is+1+qis+a > is+a+1k=qka1+ka+qk1 note that for k such that 0ka1 , the identification gives no information on the values of the k . these formulas show that all the k depend linearly on the a first states , k , k=0,  ,a1 . hence , using ( 4 ) , the generating function of l is only a function of the a first states . dividing ( 4 ) by ( 1u)(1uq ) , we can find :rar linear in the first a coordinates such that:(a.3)(1++ua1)e(ul)=(0, the a1 roots of r(u)=1++ua1 are the uk = e2ik / a for 1ka1 , which are in the unit disk . therefore , the vector ( 0,  ,a1 ) is solution of the system:(a.4)(0, if the aa system ( a.4 ) to which we added the normalization condition is invertible , then there exists a unique solution ( 0,  ,a1 ) . having the vector ( 0,  ,a1 ) , the generating function e(l ) follows from ( a.3 ) . for instance k , for k in a+1,is , is a linear combination ( independent of k ) of the k - th powers of the roots of the ( conjugate of the ) characteristic polynomial cxa+1(1+c)xa+1 , where c=2q / p(0, ) . this last property is also true for the characteristic polynomials of the relations for the k for k in is+2,is+a or > is+a+1 . as the quartic equation has explicit solutions , the expressions of the k , ka are fully explicit if a+14 + 1 , that is a4 . in particular for a=1 , the ( k)k1,is and ( k)k > is are two geometric progressions:(a.5)k1,is,k = ck0,k > is,k = p(1p)k(2p)is+10 . first , we prove that both markov chains defined by ( 2 ) , ( 3 ) are ergodic and , second , we derive formula ( 4 ) . the ergodicity is a direct consequence of foster - lyapunov criteria ( corollary 8.7 p. 214 in robert ( 2003 ) or proposition 1.3 in hairer ( 2016 ) , for instance ) . the chains are time - homogeneous and , for l0>a : e(l1l0)=e(ab0+c0g0)=a/2+f(l0)e(g0),and this last quantity is negative for large enough values of l0 because f tends to zero , either f is f(l)=(1+(lls))1 or f(l)=1lis . this proof of ergodicity works for any function f having a limit l at infinity which satisfies l < moreover , the markov chains ( ln ) defined by ( 2 ) , ( 3 ) are irreducible and aperiodic . therefore , in both cases , there exists a unique equilibrium distribution , denoted l. for the model ( 3 ) , it is characterized by the fact that l and ( lab0)++g01{lis } have the same probability generating function . to establish formula ( 4 ) we distinguish the three regimes : l<ab0 , ab0lis , l>is , and the cases b0=0 or b0=1 . for u such that |u|1 , we obtain:(a.1)2e(ul)=2e(u(lab0)++g01lis)=2e(ug01{l<ab0})+2e(ulab0+g01{ab0lis})+2e(ulab01{l>is})=e(ug0)[p(l<a)+p(l<0)]+e(ug0)[e(ula1{alis})+e(ul1{0lis})]+[e(ula1{l>is})+e(ul1{l>is})]=e(ug0)sa1(1)+e(ug0)[ua(sis(u)sa1(u))+sis(u)]+[ua(e(ul)sis(u))+(e(ul)sis(u ) ) ] . from ( a.1 ) , we get : e(ul)(1ua)=sis(u)(e(ug0)1)(1+ua)+e(ug0)(sa1(1)uasa1(u ) ) . as the probability generating function of a geometric distribution is explicitly given by e(ug0)=p(1uq)1 , with q=1p , we obtain ( 4 ) after multiplication by ( 1uq)ua . in order to compute the k for all k , we identify the coefficients of the power series of each side of ( 4 ) distinguishing cases for the values of k. for simplicity we set k=0 for k<0 . using the identities : ( ua1)(1uq)=qua+1+ua+qu1 , 
 ( u1)(ua+1)=ua+1ua+u1 , we get:(a.2)k=0(qka1+ka+qk1k)uk=k=0k(quk+a+1+uk+a+quk+1uk)=(qua+1+ua+qu1)e(ul)=q(ua+1ua+u1)sis(u)+psa1(1)uapsa1(u)=q[k = a+1is+a+1ka1ukk = ais+akauk+k=1is+1k1ukk=0iskuk]+psa1(1)uapk=0a1kuk . the following table gives the recurrence relations obtained after identification in the coefficients of ( a.2).tablekrelation0,a1aa=21ppoa1k=1ka+1,isk=21ppka1 + 2ppkais+1is+1=2qisa+(1+q)is+1ais+2,is+ak=2qka1+(1+q)ka+qk1is+a+1is+a+1=2qis+is+1+qis+a > is+a+1k=qka1+ka+qk1 note that for k such that 0ka1 , the identification gives no information on the values of the k . these formulas show that all the k depend linearly on the a first states , k , k=0,  ,a1 . hence , using ( 4 ) , the generating function of l is only a function of the a first states . dividing ( 4 ) by ( 1u)(1uq ) , we can find :rar linear in the first a coordinates such that:(a.3)(1++ua1)e(ul)=(0, the a1 roots of r(u)=1++ua1 are the uk = e2ik / a for 1ka1 , which are in the unit disk . therefore , the vector ( 0,  ,a1 ) is solution of the system:(a.4)(0, if the aa system ( a.4 ) to which we added the normalization condition is invertible , then there exists a unique solution ( 0,  ,a1 ) . having the vector ( 0,  ,a1 ) , the generating function e(l ) follows from ( a.3 ) . for instance k , for k in a+1,is , is a linear combination ( independent of k ) of the k - th powers of the roots of the ( conjugate of the ) characteristic polynomial cxa+1(1+c)xa+1 , where c=2q / p(0, ) . this last property is also true for the characteristic polynomials of the relations for the k for k in is+2,is+a or > is+a+1 . as the quartic equation has explicit solutions , the expressions of the k , ka are fully explicit if a+14 + 1 , that is a4 . in particular for a=1 , the ( k)k1,is and ( k)k > is are two geometric progressions:(a.5)k1,is,k = ck0,k > is,k = p(1p)k(2p)is+10 . to rigorously compare the variance of the simplified model ( 3 ) with the one of ( 2 ) , we choose is so that the ceiling function of the mean of the equilibrium of ( 3 ) 
 e(l) is the same as the one for ( 2 ) , i.e. 342  bp in xu et al . ( 2013):a=7,p=0.026,ls=90,=0.045,where  is a fitting parameter describing telomerase recruitment to telomeres ( xu et al . , 2013 ) ( table e.1 ) . in order to compute the corresponding mean of the equilibrium , for each choice of the aim of this section is to establish formula ( 5 ) . before entering the details just note that taking a=1 is not strictly equivalent to make a change in the length unit as the lenghts of the telomeres are integer - valued . but taking a>1 only make the results more difficult to state without exhibiting a new behavior . for the sake of simplicity , we drop the superscript i in this section and start by studying a typical couple ( xn , yn ) . the shortening of these two telomeres can be mathematically translated into the following model:(c.6)(xn+1yn+1)=((xnabn)+(yna(1bn))+),where bn is a bernoulli random variable of parameter 1/2 , and ( x0,y0)dist.. this process is an oriented simple random walk on z2 until one of the coordinates reaches zero , and can be written explicitly:(c.7)xn = x0ak=1nbk = x0ab(n,1/2),yn = y0ak=1n(1bk)=y0a(nb(n,1/2)),where b(n,1/2 ) is a binomial distribution of parameters n and 1/2 . in this case , let us define the first time one of the coordinates reaches zero , t , as : t1=inf{n0,min(xn , yn)<0}then , from ( c.7 ) and for a=1 , because ( xn ) and ( yn ) are non - increasing , we get:(c.8)p(t1>n)=p(xn0,yn0)=p(ny0b(n,1/2)x0)=k+lnk , l0(x0=k , y0=l)2nt = nlk(nt).from here , we easily derive the distribution of the time of senescence by considering all 16 independent pairs of telomeres:(c.9)p(t > n)=p({senescenceatthenthgeneration})=p(min1k32lkn0)=p(i1,16,min(xni , yni)0)=p(min(xn , yn)0)16=p(t1>n)16 . formulas ( c.8 ) , ( c.9 ) lead to ( 5 ) , which gives the expected time of senescence ( 6 ) using e(t)=n=0p(t > n ) . the aim of this section is to prove formulas ( 9 ) , ( 10 ) . to do so we first note that , using a generalized version of the central limit theorem , for any fixed positive number t , and any sequence untn:(d.10)b(un,12)12unnn+l12n(0,t),where b(un,12 ) is a binomial distribution with parameters un and 1/2 . as in ( c.8 ) , we get , with n(x0)=2x0wx02x0:p(2x0tx01x0<w)=p(tx01>2x0wx0)=p(tx01>n(x0))=p(n(x0)x0b(n(x0),12)x0)=p(|b(n(x0),12)n(x0)2x0|x0n(x0)2x0 ) . as x0 tends to  , thanks to ( d.10 ) , we obtain : p(2x0tx01x0w)p(|n|2w2),where nn(0,2 ) . the speed of convergence is given by the berry - essen theorem in terms of the cumulative distribution function of 2x0tx01x0 , say fx0:xr,x0n|fx0(x)erfx|cx01/2,where c is an absolute constant . see table e1parameters used in this study.table e1alength of the 3-end overhang.sthreshold length of the shortest telomere inducing senescence.lslength threshold for the function f , below which telomerase is recruited to the telomere with probability 1.islength threshold for the simplified model of telomerase recruitment.pparameter of the geometric random variable gn.parameter of the function f , fitted on experimental data ( xu et al . , 2013 ) .
replicative senescence , induced by telomere shortening , exhibits considerable asynchrony and heterogeneity , the origins of which remain unclear . here , we formally study how telomere shortening mechanisms impact on senescence kinetics and define two regimes of senescence , depending on the initial telomere length variance . 
 we provide analytical solutions to the model , highlighting a non - linear relationship between senescence onset and initial telomere length distribution . 
 this study reveals the complexity of the collective behavior of telomeres as they shorten , leading to senescence heterogeneity .