Tumor growth, $R$-positivity, Multitype branching and Quasistationarity

TL;DR

This paper investigates the mathematical properties of tumor growth models using Markov processes and positive matrices, establishing conditions for their long-term behavior and applying these to understand tumor size dynamics.

Contribution

It introduces new conditions for R-positivity in Markov processes and positive matrices, and applies these to analyze tumor size asymptotics in supercritical regimes.

Findings

Conditions for R-positivity established

Asymptotic behavior of tumor sizes derived

Application to supercritical tumor growth models

Abstract

Motivated by tumor growth models we establish conditions for the $R-$positivity of Markov processes and positive matrices. We then apply them to obtain the asymptotic behavior of the tumors sizes in the supercritical regime.

Equations

$$\frac{d{x}_p}{dt}=h(x_p),x_p(0)=1.$$

$$\frac{\partial u}{\partial t}(x,t)+\frac{\partial (hu)}{\partial x}(x,t)=0,x>1,t>0,$$

$$h(1)u(1,t)={\int }_1^{\infty }\beta (x)u(x,t)dx+\beta (x_p(t)).$$

$$\mathcal{L}f(\eta )=\sum _{x\in \Lambda }\sum _{y\in {\Lambda }_0}q(x,y)\eta (x)[f({\eta }^{x,y})-f(\eta )]+\sum _{x\in \Lambda }\beta (x)\eta (x)[f({\eta }^{+})-f(\eta )].$$

$${\eta }^{x,y}(z)=\{\begin{array}{ll}\eta (z),& z\notin \{x,y\},\\ \eta (x)-1,& z=x\\ \eta (y)+1& z=y,y\ne 0.\end{array}{\eta }^{+}(z)=\{\begin{array}{ll}\eta (z),& z\ne 1,\\ \eta (1)+1,& z=1.\end{array}$$

$${\kappa }_0:={\int }_0^{\infty }{\gamma }_1(t)dt>1.$$

$${\gamma }_1(t)=\sum _x\beta (x){\mathbb{E}}_1(1\{Y_t=x\})$$

$$\underset{t\to \infty }{\lim }\frac{{\eta }_t}{|{\eta }_t|}=\nu ,\text{in probability}.$$

$$\dot{V}(x)=\sum _{y\in {\Lambda }_0}a(x,y)V(y).$$

$$\frac{d}{dt}{g}_x(t)⩽{\mathbb{E}}_x\dot{V}(X_t)1\{X_t\notin {\mathcal{U}}_1\}.$$

$$\frac{d}{dt}{g}_x(t)⩽(\log \rho ){\mathbb{E}}_{x}V(X_t)1\{X_t\notin {\mathcal{U}}_1\}=(\log \rho )g_x(t).$$

$${\mathbb{P}}_x$$

$$=\frac{V(x)}{\tilde{K}}{\mathbb{E}}_x\left[\frac{V(X_1)}{V(x)}\cdots \frac{V(X_{n-1})}{V(X_{n-2})}\frac{V(X_n)}{V(X_{n-1})}1\{X_1\in {\mathcal{U}}_1^c,\dots ,X_n\in {\mathcal{U}}_1^c\}\right]$$

$$⩽{\mathbb{E}}_x\left[\frac{V(X_1)}{V(x)}\cdots \frac{V(X_{n-1})}{V(X_{n-2})}\frac{V(X_n)}{V(X_{n-1})}1\{X_1\in {\mathcal{U}}_1^c,\dots ,X_n\in {\mathcal{U}}_1^c\}\right]$$

$$\displaystyle=\mathbb{E}_{x}\left\{\mathbb{E}\left[\frac{V(X_{1})}{V(x)}\cdots\frac{V(X_{n-1})}{V(X_{n-2})}\frac{V(X_{n})}{V(X_{n-1})}{{\bf\rm 1}}\{X_{1}\in\mathcal{U}_{1}^{c},\dots,X_{n}\in\mathcal{U}_{1}^{c}\}\ \Bigg{|}\ X_{1},\dots,X_{n-1}\right]\right\}$$

$$\displaystyle\leqslant\mathbb{E}_{x}\left[\frac{V(X_{1})}{V(x)}\cdots\frac{V(X_{n-1})}{V(X_{n-2})}{{\bf\rm 1}}\{X_{1}\in\mathcal{U}_{1}^{c},\dots,X_{n}\in\mathcal{U}_{1}^{c}\}\right]\cdot{\sup_{x^{\prime}\in\mathcal{U}_{1}^{c}}\mathbb{E}\big{(}\textstyle\frac{V(X_{n})}{V(X_{n-1})}\big{|}X_{n-1}=x^{\prime}\big{)}}$$

$$⩽\rho \cdot {\mathbb{E}}_x\left[\frac{V(X_1)}{V(x)}\cdots \frac{V(X_{n-1})}{V(X_{n-2})}1\{X_1\in {\mathcal{U}}_1^c,\dots ,X_n\in {\mathcal{U}}_1^c\}\right]$$

$$⩽\cdots ⩽{\rho }^{n-1}\cdot {\mathbb{E}}_x\left[\frac{V(X_1)}{V(x)}\right]⩽{\rho }^{n-1}\cdot \underset{x\in \Lambda }{\sup }{\mathbb{E}}_x\left[\frac{V(X_1)}{V(x)}\right]⩽C{\rho }^{n-1}.$$

$$a(z)=\underset{x\in \Lambda \setminus \{z\}}{\inf }a(x,z),\alpha =\sum _{z\in \Lambda }a(z),C=\underset{x\in \Lambda }{\sup }a(x,0).$$

$$\underset{t\to \infty }{\lim }-t^{-1}\log {\mathbb{P}}_x(X_t=y)=\gamma ,$$

$$-t^{-1}\log {\mathbb{P}}_x(X_t=x)⩾\gamma .$$

$$\gamma ⩽-t^{-1}\log {\mathbb{P}}_x(X_t=x)⩽-t^{-1}\log {\mathbb{P}}_x(X_t\ne 0)⩽C.$$

$$\rho =e^{-\alpha }<e^{-\gamma }=R.$$

$${\mathbb{P}}_x({\tau }_{{\mathcal{U}}_1}>n,{\tau }_0>n)⩽{\mathbb{P}}_x({\tau }_{{\mathcal{U}}_1}>n)⩽e^{-\alpha n}={\rho }^n.$$

$$\underset{n\to \infty }{\lim }{\mathbb{P}}_x(X_n=y|{\tau }_0>n)=\nu (y),\underset{n\to \infty }{\lim }{R}^n{\mathbb{P}}_x({\tau }_0>n)=\mu (x).$$

$${\mathbb{P}}_x({\tau }_0>n)$$

$$⩽{\mathbb{P}}_x({\tau }_{{\mathcal{U}}_1}>n)+{\mathbb{E}}_x({\mathbb{P}}_{X_{{\tau }_{{\mathcal{U}}_1}}}({\tau }_0>n-{\tau }_{{\mathcal{U}}_1})1\{{\tau }_{{\mathcal{U}}_1}⩽n\})$$

$$⩽e^{-\alpha n}+{\mathbb{E}}_x(g(n-{\tau }_{{\mathcal{U}}_1})),$$

$$\mu (x)=\underset{n\to \infty }{\lim }{R}^{-n}{\mathbb{P}}_x({\tau }_0>n)$$

$$⩽\underset{n}{\mathrm{limsup}}{\mathbb{E}}_x\left(R^{-(n-{\tau }_{{\mathcal{U}}_1})}g(n-{\tau }_{{\mathcal{U}}_1})R^{-{\tau }_{{\mathcal{U}}_1}}\right)$$

$$⩽\left(\underset{n}{\sup }{R}^{-n}g(n)\right){\mathbb{E}}_x{R}^{-{\tau }_{{\mathcal{U}}_1}}=:C_1.$$

$${\mathbb{E}}_x(R^{-{\tau }_{{\mathcal{U}}_1}})=\sum _{n\in \mathbb{N}}{R}^{-n}\mathbb{P}({\tau }_{{\mathcal{U}}_1}=n)⩽\sum _{n\in \mathbb{N}}{R}^{-n}\mathbb{P}({\tau }_{{\mathcal{U}}_1}>n-1)=\frac{1}{R}\sum _{n\in {\mathbb{N}}_0}{R}^{-n}{\rho }^n<\infty .$$

$$\mu (x)$$

$$⩾\gamma \alpha \underset{y\in {\mathcal{U}}_1}{\min }\mu (y)=:c_1>0.$$

$$\sum _{y\in \Lambda }{\mathbb{E}}_y((\sum _{x\in \Lambda }{Z}_1(x)\mu (x){)}^2)\nu (y)<\infty .$$

$$\underset{n\to \infty }{\lim }{\mathbb{E}}_x{\left[R^{-n}\sum _{y\in \Lambda }{Z}_n(y)f(y)-W\sum _{y\in \Lambda }\nu (y)f(y)\right]}^2=0.$$

Full text

1. Introduction

Our work in this paper explores the $R-$positivity and quasistationarity for a class of models motivated by the tumor growth model studied by L. Triolo in [Tri05]. His work in turn was motivated by works of [IKS00, Str03] where a continuum model of tumor growth is studied. In that model a primary tumor starts somewhere (space is not a variable) with one cell and grows at a rate $h(x)$ where $h$ is taken to be a Gompertzian law. That is, the number of cells of the primary tumor $x_{p}$ is given by the following ordinary differential equation

[TABLE]

Here $h(x_{p})=ax_{p}\log(\frac{N}{x_{p}})$, thus asymptotically the cell size approaches $N$. The malignant behavior is modeled by creation of one new (metastases) cell at a rate $\beta(x)$ which is an increasing function of $x$. Each new cell grows and proliferates according to the same rule as the primary tumor. The proliferation rate $\beta(x)$ is taken to be $\beta(x)=\kappa x^{r}$ where $0<r\leqslant 1$ and $\kappa$ is a constant called the colonization constant. The evolution of the distribution of metastases $u(x,t)$ is given by the following equation.

[TABLE]

with a boundary condition at $x=1$ given by

[TABLE]

This equation is analyzed in [IKS00, Str03], where it is shown that $u$ increases asymptotically exponentially in $t$.

In [Tri05] a birth and death model was proposed with a view towards including the random localized (proliferation occuring at size one) nature of cell growth as well as to include the possibility of modeling immune response by a nonzero probability of death when the cell size is one. It is a microscopic model in which individual particles (tumor sizes) evolve independently as a birth and death process on $\mathbb{N}\cup\{0\}$ with [math] being an absorbing state. The birth and death rates are chosen so that the drift matches the Gompertzian law of the continuum model. If we denote the occupation number of site $x$ by $\eta(x)$, then with rate $\sum_{1}^{\infty}\beta(x)\eta(x)$ a particle is created at site one (metastases). For such a model it is shown in [Tri05] that as $t\to\infty$ the expected occupation numbers converge to 1) zero 2) a constant nonzero value or 3) diverge exponentially to infinity depending on the value of a parameter $\kappa_{0}$ being 1) less than one 2) equal to one 3) greater than one. The parameter $\kappa_{0}$ is the expected total number particles created by one particle before being eliminated and is defined precisely below.

We study the asymptotic behavior of this model in the supercritical regime $\kappa_{0}>1$. The main tool is to identify the process with a multitype branching process and to show that its mean matrix is $R-$positive, which leads to a Kesten-Stigum type theorem for the asymptotic behavior of the distribution of tumor sizes [Moy67, Eng07, Eng15]. In the course of the proof we provide a general criteria to establish the $R-$ positivity of positive matrices by identifying certain transformations of them with rates matrices of absorbed Markov processes.

1.1. The Model.

We use an interacting particle system to model our process where the particles move according to a Markov process with a countable state space $\Lambda$ and get absorbed at the state [math]. We denote $(\eta_{t}(x),\,x\in\Lambda)$ the number of particles at site $x$ at time $t$.

Let $Q_{0}=(q(x,y),\,x,y\in\Lambda_{0})$ be a rates matrix of a pure jump Markov process on $\Lambda_{0}:=\Lambda\cup\{0\}$. We assume that $Q$ is irreducible, [math] is absorbing ($q(0,y)=0$ for all $y\in\Lambda$) and denote $Q$ the restriction of $Q_{0}$ to $\Lambda$. We use the convention $q(x,x)=-q(x)=-\sum_{y\in\Lambda_{0}\setminus\{x\}}q(x,y)$. Observe that $q(x,0)=-\sum_{y\in\Lambda}q(x,y)$ is the absorption rate from state $x$. The process can be described as follows.

1. (1)

Each particle evolves independently according to the rates matrix $Q=(q(x,y),\,x,y\in\Lambda)$. In particular, a particle a site $x$ is absorbed at a state that we call 0 at rate $q(x,0)$. 2. (2)

Particles are created at a state that we call 1 at rate $\sum_{x}\beta(x)\eta_{t}(x)$.

More precisely, $(\eta_{t},t\geqslant 0)$ is a Markov process with state space $\mathbb{N}^{\Lambda}$ and generator given by

[TABLE]

Here

[TABLE]

Of particular interest is the case where $Q$ is a birth and death process on $\mathbb{N}$ with absorption at [math] and birth and death rates given by $\lambda_{x}=ax\ln(N+1)$ and $\mu_{x}=ax\ln(x+1)$ respectively. We denote this process by $(Y_{t})_{t\geqslant 0}$. With such a choice the drift $\lambda_{x}-\mu_{x}$ follows a Gompertzian law $g(x)=ax\ln\frac{N+1}{x+1}$. Here $x$ represents the tumor size. Also of special interest is the case $\beta(x)=\kappa x^{r}$ for $0<r\leqslant 1$. In [Tri05] it is proved that in this situation the process is supercritical if and only if

[TABLE]

Here

[TABLE]

is the expected creation rate of a single particle. More precisely, it is proved that

1. (a)

$\kappa_{0}<1\Rightarrow\mathbb{E}\eta_{t}(x)\to 0$ for all $x\in\mathbb{N}$. 2. (b)

$\kappa_{0}=1\Rightarrow\mathbb{E}\eta_{t}(x)$ converges to a constant. 3. (c)

$\kappa_{0}>1\Rightarrow\mathbb{E}\eta_{t}(x)$ grows exponentially fast.

The proof can be extended with no difficulty to our more general situation. It can also be obtained from the following observation.

We can couple the total number of particles with a single-type Galton-Watson process. For a given particle, consider the total number of children given birth by a particle during its lifetime. This is the offspring distribution. In this way the total number of particles of generation $n$ is a Galton-Watson process (but not the total number of particles at time $t$). From this observation, and by means of irreducibility, (a), (b), (c) can be extended to

1. (d)

$\kappa_{0}\leqslant 1\Rightarrow\eta_{t}(x)\to 0$ for all $x\in\mathbb{N}$. 2. (e)

$\kappa_{0}>1\Rightarrow\eta_{t}(x)$ grows exponentially fast for all $x\in\mathbb{N}$.

Both statements hold almost surely.

The process can be constructed in a standard way and so, we omit the proof of its existence.

After proving general conditions to guarantee $R-$positivity of positive matrices and the asymptotic behavior of supercritical multitype branching process we will obtain for this model the following.

Theorem 1.2**.**

Under adequate assumptions on $\beta$ and $Q$, if $\kappa_{0}>1$ there is positive probability of non-extinction and on this event we have

[TABLE]

Here $\nu$ is a finite measure, which is the left eigenvector of a matrix obtained as a transformation of $Q$ (see section 3 for details). The paper is organized as follows. In Section 2 we establish two general criteria to prove $R-$positivity of nonnegative matrices. In both of them we transform the matrix to obtain a sub-Markovian operator. Next we apply a Lyapunov type criteria in the first case (Theorem 2.1) and a Doëblin type argument in the second one, Theorem 2.5. In Section 3 we apply these criteria to prove Theorem 1.2.

Similar strategies have been applied in a series of papers by N. Champagnat, D. Villemonais ([CV16, CV17, CV19] among them) to establish existence of quasi-stationary distributions and uniform convergence towards them. Our techniques differ from theirs. They are probabilistic in nature and based on a theorem by Kesten, Ferrari and Martínez [FKM96] that provides probabilistic conditions to establish $R-$positivity.

2. Lyapunov functions, Döeblin conditions and $R$-positivity

In this section we consider continuous or discrete time Markov chains and give conditions under which the process is $R$-positive. We will use them in Section 3 to establish the asymptotic behavior of the tumor growth model. The processes in this section should not be confused with the driving process with rates $Q$ of the previous section. For that reason, in this section we use the letters $A$, $A_{0}$, $\tilde{A}$ for generators instead of $Q$.

Let $X=(X_{t},t\geqslant 0)$ be a continuous time pure jump Markov process with rates matrix $A_{0}=(a(x,y),x,y\in\Lambda\cup\{0\})$. We use again the convention $a(x,x)=-\sum_{y\neq x}a(x,y)$. We assume that $X$ is absorbed at zero and denote $A$ the restriction of $A_{0}$ to $\Lambda$, so that $-\sum_{y\in\Lambda}a(x,y)=a(x,0)$ is the absorption rate from state $x$. For a function $V\colon\Lambda\to\mathbb{R}_{\geqslant 0}$, define the drift of $V$ at $x\in\Lambda$ by

[TABLE]

We always assume $V(0)=0$ and that the sum above is well defined and finite.

Theorem 2.1** (Lyapunov condition).**

Assume there exists $V\colon\Lambda\to\mathbb{R}_{\geqslant 0}$ with $V(x)\to\infty$ as $x\to\infty$ and such that $\lim_{x\to\infty}\frac{\dot{V}(x)}{V(x)}=-\infty$. Then $X$ is $R$-positive for some $0<R\leqslant 1$.

Remark 2.2*.*

This theorem has to be compared with [CV17, Theorem 5.1]. Our assumptions seem to be more restrictive than those on [CV17], however our proof can be easily modified to fit those assumptions. We prefer to write it in this way since this is a condition that can actually be checked in practice. Also, although not exactly the same, the conclusions are similar (but [CV17] is more general and includes and exponential convergence to quasistationarity statement). We decided to include our proof since it is short and simple. It is based on ideas and the following theorem from [FKM96].

Theorem 2.3** ([FKM96, Theorem 1]).**

Assume $P=\exp(A)$ is irreducible and call $R$ its decay parameter. Let $\tau=\inf\{n\in\mathbb{N}\colon X_{n}=0\}$. Suppose that there exist a finite set $\mathscr{U}_{1}\subset\Lambda$, a state $x^{\prime}\in\mathscr{U}_{1}$, $\rho<R$, and a positive constant $\kappa$ such that for all $x\in\mathscr{U}_{1}$ and $n\geqslant 0$, $\mathbb{P}_{x}(\tau>n;\,X_{1}\notin\mathscr{U}_{1},\dots,X_{n}\notin\mathscr{U}_{1})\leqslant\kappa\rho^{n}$; Then $P$ is $R$-positive and its left eigenvector $\nu$ is summable.

Remark 2.4*.*

The statement of [FKM96, Theorem 1] is more involved since their subset $\mathscr{U}_{1}$ can be infinite. In our case, conditions (1.23) and (1.24) in [FKM96] follow immediately from irreducibility of $A$. The hypotheses of the theorem above correspond to what they call (1.22). In our statement $R$ is what they call $1/R$. The main ingredient of the proof of Theorem 2.1 is to show that this condition is in fact verified.

Proof of Theorem 2.1.

Let $X=(X_{t},t\geqslant 0)$ be a process with a generator $A$. Given $x\in\Lambda$ define $g_{x}(t)=\mathbb{E}_{x}(V(X_{t}){{\bf\rm 1}}\{X_{t}\notin\mathscr{U}_{1}\})$. We have for $x\notin\mathscr{U}_{1}$

[TABLE]

Given $0<\rho<1$, we can choose $\mathscr{U}_{1}$ such that $-\dot{V}(x)/V(x)>-\log\rho$ if $x\notin\mathscr{U}_{1}$. Then, for such $\mathscr{U}_{1}$ we have

[TABLE]

This and $g_{x}(0)=V(x)$ gives us $g_{x}(t)\leqslant V(x)\rho^{t}$ and in particular $\mathbb{E}_{x}(V(X_{1}))/V(x)\leqslant\rho$ for any $x\notin\mathscr{U}_{1}$ and also $\sup_{x}\mathbb{E}_{x}(V(X_{1}))/V(x)\leqslant C$ for every $x\in\Lambda$. Given $\tilde{K}$, we can enlarge $\mathscr{U}_{1}$ if necessary to ensure that $x\notin\mathscr{U}_{1}$ implies $V(x)\geqslant\tilde{K}$. Then we have for any $x$ with $V(x)\leqslant K$,

[TABLE]

By means of Theorem 2.3, we get that discrete time chain $(X_{n})$ is $R-$ positive, which is equivalent to the $R-$positivity of its transition probability matrix ${\rm exp}(A)$. Finally, observe that there is a bijection between the eigenvalues of $A$ and ${\rm exp}(A)$ and that the eigenvectors are exactly the same. This gives us the $R-$positivity of $A$, or equivalently, the process $(X_{t})_{t\geqslant 0}$. ∎

Next we give a different condition that guarantees $R$-positivity. It also has to be compared with [CV17, Remark 11] and [FM07, Theorem 1.4]. Define

[TABLE]

Theorem 2.5**.**

If $\alpha>C$, then $A$ is $R$-positive. Moreover, the right eigenvector $\mu$ verifies $0<c_{1}\leqslant\mu(x)\leqslant C_{1}$ for some positive constants $c_{1}$ and $C_{1}$ and all $x\in\Lambda$. As a consequence, $\nu$ is summable.

Remark 2.6*.*

Under this condition, Ferrari and Marić [FM07] proved existence of a quasi-stationary distribution (QSD) and Jacka and Roberts uniqueness [JR95]. Here we prove $R$-positivity, which also leads to existence of such a QSD.

Proof.

The argument is reminiscent to the one in [CV16, Theorem 4.1], although applied in a different situation. Let $(X_{t},t\geqslant 0)$ be a Markov process with transition rates $A$. By [Kin63b, Theorem 1], the parameter $\gamma$ defined by

[TABLE]

is well defined, independent of $x,y$ and verifies for every $x\in\Lambda,\,t\geqslant 0$,

[TABLE]

Observe that for every $x\in\Lambda$ we have $\mathbb{P}_{x}(X_{t}\neq 0)\geqslant e^{-Ct}$. Hence

[TABLE]

Without loss of generality, we can assume that there exists a finite set $\mathscr{U}_{1}\subset\Lambda$ such that $\sum_{y\in\mathscr{U}_{1}}\inf_{x\in\Lambda}a(x,y)=\alpha>C.$ Then, for $R=e^{-\gamma}$, we can take $\rho>e^{-\alpha}$ such that

[TABLE]

Let us check that $P:=\exp(A)$ verifies the hypothesis of Theorem 2.3. Observe that since $\mathbb{P}_{x}(\tau_{\mathscr{U}_{1}}>n)=\mathbb{P}_{x}(X_{k}\notin\mathscr{U}_{1},\mbox{for }k\in\{1,2,\ldots,n\})\leqslant e^{-\alpha n}$, we have

[TABLE]

Then, $P$ is $R$-positive and there exist left and right eigenvectors $\nu$, $\mu$ respectively with eigenvalue $1/R$. By [JR95], $\nu$ is the only QSD of $P$. Also, $\sum_{x\in\Lambda}\nu(x)<\infty$ and by [SVJ66, Theorem 3.1] the following limits hold

[TABLE]

Next, we compute the bounds for $\mu$. Notice that

[TABLE]

where $g(k)=\max_{y\in\mathscr{U}_{1}}\mathbb{P}_{y}(\tau_{0}>k)$. Then,

[TABLE]

The first factor is finite since $R^{-n}g(n)$ has a limit which does not depend on $x$. To bound the second factor, observe that since $\mathbb{P}_{x}(\tau_{\mathscr{U}_{1}}>n)\leqslant e^{-\alpha n}=\rho^{n}$ for all $n\geqslant 1$, we have

[TABLE]

For the lower bound we compute,

[TABLE]

∎

3. Tumor growth and multitype branching

In this section we will address the asymptotic behavior of $(\eta_{t})$ as time goes to infinity in the supercritical regime. Our proofs will strongly rely on a theorem by Moy (a version of Kesten-Stigum theorem for countable types) for the behavior of supercritical multitype branching processes. We state the theorem for completeness. We will slightly abuse notation by using $\mathbb{E}_{y}$ and $\mathbb{P}_{y}$ to denote that the process $Z$ starts with one particle at $y$ at $n=0$. We use $|Z_{n}|$ for the total number of particles in the sytem at time $n$, that is $|Z_{n}|=\sum_{y\in\Lambda}Z_{n}(y)$.

Theorem 3.1** ([Moy67, Theorem 1]).**

Let $(Z_{n})_{n\geqslant 0}$ be a multitype branching process with a countable types space $\Lambda$ and mean matrix $M$ with entries $m(x,y)=\mathbb{E}_{x}(Z_{1}(y))$, $x,y\in\Lambda$. Assume $M$ is aperiodic, irreducible and $R$-positive with normalized left and right eigenvectors $\nu,\,\mu$ respectively, $R>1$ (supercritical) and

[TABLE]

Then there is a real valued random variable $W$, with $\mathbb{E}(W^{2})<\infty$ such that for every $f$ with $\sup_{x}f(x)/\mu(x)<\infty$ we have

[TABLE]

In particular, for every $y\in\Lambda$, $\mathbb{E}_{x}[R^{-n}Z_{n}(y)-W\nu(y)]^{2}\to 0$ and if $\inf_{x\in\Lambda}\mu(x)>0$, $\mathbb{E}_{x}[R^{-n}|Z_{n}|-W]^{2}\to 0$.

As a consequence, we obtain the following corollary, as in [JS17]

Corollary 3.3**.**

In the conditions of Theorem 3.1, assume that for every $x\in\Lambda$ and $Z_{0}=\delta(x,\cdot)$, conditioned on non-extinction, there is $y\in\Lambda$ with $Z_{n}(y)\to\infty$ in probability as $n\to\infty$, then $Z_{n}/|Z_{n}|\to\nu$ in probability on this event.

Before embarking on the proof, we need to introduce some notation. Define the moment generating function $\mathbf{f}\colon[0,1]^{\Lambda}\to\mathbb{R}^{\Lambda}$ by

[TABLE]

If we define

[TABLE]

we get, as for single-type Galton-Watson processes, $\mathbf{f}_{n+1}(\mathbf{s})=\mathbf{f}(\mathbf{f}_{n}(\mathbf{s}))$. We also define $\mathbf{q}:=(q_{x},\,x\in\Lambda)$, with $q_{x}$ being the absorption probability when the process starts with one individual of type $x$, i.e. $q_{x}:=\mathbb{P}_{x}\left(\bigcup_{n\in\mathbb{N}}\{Z_{n}=\mathbf{0}\}\right)$.

Proof.

The proof follows mainly [JS17]. The function $\mathbf{f}$ has at least two fixed points: $\mathbf{1}=(1,x\in\Lambda)$ and $\mathbf{q}$. It is also known that if $\mathbf{s}$ is a fixed point of $\mathbf{f}$ different from $\mathbf{1}$, then we have $s_{x}<1$ for every $x\in\Lambda$.

For the first claim, $\mathbf{1}$ is clearly a solution to $\mathbf{f}(\mathbf{s})=\mathbf{s}$ and for the extinction event, we have

[TABLE]

That is, $\mathbf{q}=\displaystyle\lim_{n\to\infty}\mathbf{f}_{n}(\mathbf{0})$. Since $\mathbf{f}_{n+1}(\mathbf{0})=\mathbf{f}(\mathbf{f}_{n}(\mathbf{0}))$ and $\mathbf{f}$ is continuous in this topology, we get $\mathbf{q}=\mathbf{f}(\mathbf{q})$.

For the second claim, let $\mathbf{s}\in[0,1]^{\Lambda}$ be a fixed point of $\mathbf{f}$ with $s_{x}<1$ for some $x$ and assume there is $y\in\Lambda$ with $f_{y}(\mathbf{s})=1=s_{y}$. Then $\mathbf{s}$ is a fixed point of $\mathbf{f}_{n}$ for every $n$. Let $\mathbf{t}\in[0,1]^{\Lambda}$ be such that

[TABLE]

[TABLE]

Since $\mathbf{s}\leqslant\mathbf{t}$, $f_{ny}(\mathbf{s})=1$, and the function $f_{ny}$ is monotone increasing, we get $1=f_{ny}(\mathbf{t})$ and hence

[TABLE]

Then $\mathbb{P}_{y}(Z_{n}=\mathbf{z})=0$ if $z_{x}>0$, since $\sum_{\mathbf{z}\in\mathbb{N}^{\Lambda}}\mathbb{P}_{y}(Z_{n}=\mathbf{z})=\mathbf{1}.$ So,

[TABLE]

for every $n\in\mathbb{N}$. This contradicts the irreducibility of $M$. Then $\mathbf{s}<\mathbf{1}$.

Let $\mathbf{s}$ be a fixed point of $\mathbf{f}$ different from $\mathbf{1}$. We have $s_{x}<1$ for every $x\in\Lambda$. Denote $\Omega_{\rm surv}=\{|Z_{n}|\to\infty\}=\{|Z_{n}|\neq 0\,\text{ for very }n\geqslant 0\}$. Let $y$ be such that $Z_{n}(y)\to\infty$ in probability on $\Omega_{\rm surv}$, then

[TABLE]

By dominated convergence we get

[TABLE]

and

[TABLE]

and hence $s_{x}=q_{x}$ for every $x\in\Lambda$. We conclude that the only fixed points are $\mathbf{q}$ and $\mathbf{1}.$

Let $\mathbf{r}=(r_{x},\,\,x\in\Lambda)$ be the vector with coordinates $r_{x}=\mathbb{P}_{x}(W=0)$. It turns out that $\mathbf{r}$ is a fixed point of $\mathbf{f}$. To see that, we compute for $x\in\Lambda$,

[TABLE]

As a consequence, we get the dichotomy $\mathbf{r}=\mathbf{q}$ or $\mathbf{r}=\mathbf{1}$. Now, by Theorem 3.1, Taking $f=\delta_{y}$, we get

[TABLE]

and taking $f\equiv 1$ and using $\inf\mu(y)>0$, we get

[TABLE]

Then $\lim_{n\to\infty}R^{-n}|Z_{n}|=W$ in probability and hence $\{|Z_{n}|\to 0\}\subseteq\{W=0\}$. Since both events have the same probability we get

[TABLE]

This fact allow us to compute (on $\Omega_{\rm surv}$) the following limit in probability,

[TABLE]

∎

3.1. Identification with a multitype branching process

The process $(\eta_{t})$ defined by (1.1) can be identified with a multitype branching process in the following way. Recall that $q(x,x)=-\sum_{y\neq x,y\neq 0}q(x,y)$ and $q(x)=-q(x,x)$. An individual of type $x\in\Lambda$ gives birth (and die) at an exponential time of parameter $a(x):=\beta(x)+q(x)$ and the offspring distribution is given by

• •

with probability $\frac{\beta(x)}{a(x)}$, one child of type 1 and one child of type $x$.

• •

with probability $\frac{q(x,y)}{a(x)}$, one child of type $y$.

According to this, the mean matrix $M=(m(x,y),\,x,y\in\Lambda)$ for the skeleton chain is given by

[TABLE]

The expected number of individuals for the continuous time process at time $t$ is given by the matrix $\exp(tA)$, where $A$ has entries

[TABLE]

Here $\delta(x,y)=1$ if $x=y$ and [math] otherwise.

Proposition 3.4**.**

Assume $\beta$ is bounded above and one of the following conditions is verified.

1. (1)

There is $V\colon\Lambda\to\mathbb{R}_{+}$ such that

[TABLE] 2. (2)

$\sup_{x\in\Lambda}\beta(x)-\inf_{x\in\Lambda}(\beta(x)-q(x,0))<\inf_{x\in\Lambda}\beta(x)$.

Then, $A$ is $R$-positive with right eigenvector $\mu$ and left eigenvector $\nu$ such that $\sum_{x\in\Lambda}\mu(x)\nu(x)<\infty$. Moreover, $\sum_{x\in\Lambda}\nu(x)<\infty$.

As a consequence, we obtain the following corollaries.

Corollary 3.5**.**

In the notation of Proposition 3.4, if condition (2) is verified, $Q$ has a Yaglom limit (i.e. $\exp(tQ)(x,y)/(1-\exp(tQ)(x,0))\to\nu(y)$ as $t\to\infty$ for every $x\in\Lambda$) and $\kappa_{0}>1$, then (1.3) holds.

Corollary 3.6**.**

In the notation of Proposition 3.4, if $Q$ is the rates matrix of a one dimensional birth and death process that verifies,

1. (a)

$(q(x,x-1))_{x\in\mathbb{N}}$* is monotone increasing,* 2. (b)

$\lim_{x\to\infty}\frac{q(x,x+1)}{q(x,x-1)}=\ell<1,$** 3. (c)

$\sum_{x\in\mathbb{N}}\frac{1}{q(x,x-1)}<\infty.$**

and $\kappa_{0}>1$, then (1.3) holds.

Proof of Proposition 3.4.

Let $\bar{\beta}=\sup_{x}\beta(x)$ and consider the matrix $\tilde{A}$ with coefficients $\tilde{a}(x,y)=a(x,y)-\bar{\beta}\delta(x,y)$. Since $\sum_{y\neq x}\tilde{a}(x,y)\leqslant 0$ for every $x\in\mathbb{N}$, we can think of $\tilde{A}$ as the rates matrix of a process absorbed at zero that we call $X=(X_{t},t\geqslant 0)$. Observe that the absorption rate from state $x$ is given by $-\sum_{y\neq x}\tilde{a}(x,y)$. We are going to prove that $X$ is $R$-positive and hence the same holds for $\exp(tA)$ for every $t>0$. Assume condition 1. is verified, then we compute for this chain, for $x\neq 0,1$,

[TABLE]

Hence, we have $\dot{V}(x)/V(x)\to-\infty$ and $V(x)\to\infty$ as $x\to\infty$ and we can apply Theorem 2.1 to get the $R$-positivity of $\tilde{A}$, $A$ and ${\rm exp}(tA)$ for every $t>0$.

If condition 2 is verified, we apply Theorem 2.5 instead. We consider again the Markov process with rates matrix $\tilde{A}$. For this matrix we have $\alpha(1)=\inf_{x}\tilde{a}(x,1)\geqslant\inf_{x}{\beta(x)}$ and

[TABLE]

Then, we can apply Theorem 2.5.

The limit (1.3) will be a consequence of Theorem 3.1 and some additional considerations. Once $R-$ positivity is proved, we need to check (3.2). This is a delicate condition which is not simple to prove in general. We will prove that in fact (3.2) is verified under our hypotheses. Next, we also need to show that Corollary 3.3 can be applied to get (1.3).

Proof of Corollary 3.5.

Under condition (2) we have shown in Theorem 2.5 that $\mu$ is bounded above and then (3.2) reduces to

[TABLE]

Observe that the total number of particles in the system at time $t$ is stochastically dominated (uniformly in $y$) by a continuous time Galton-Watson process with binary branching and reproduction rate $\sup_{x}\beta(x)$, which has a second moment uniformly bounded in $y$. Since $\sum_{x}\nu(x)<\infty$, (3.2) holds in this case and we get for $y\in\Lambda$,

[TABLE]

Since $\inf_{x}\mu(x)>0$, we also get

[TABLE]

Now we prove that for every $x\in\Lambda$ there exists a $y\in\Lambda$ such that $P_{x}(Z_{n}(y)\to\infty)=1$. Observe that, with the exception of the ancestral particle, every other particle is born at one and hence, if it is still in the system (has not been absorbed yet), the probability of being at one is (for some $t$) $\exp(tQ)(1,1)/(1-\exp(tQ)(1,0))$, the conditional probability of being at one given that it was not absorbed, which has a positive limit and is positive for every finite time due to the existence of Yaglom limit for $Q$. Hence its infimum $c$ is larger than zero. We have for every $K$

[TABLE]

By bounded convergence theorem we have

[TABLE]

and using Borel-Cantelli’s Lemma we get $Z_{n}(1)\to\infty$ a.s. in $\Omega_{\rm surv}$ and by Corollary 3.3

[TABLE]

To go from discrete to continuous time we use [Kin63a, Theorem 2]. The following computations are conditioned on $\Omega_{\rm surv}$. Given $\epsilon>0,$ for every $y\in\Lambda$ we consider the function $g_{y}\colon[0,\infty)\to[0,1]$ defined by

[TABLE]

We show that $g_{y}$ is continuous. Let $0\leqslant s<t$. We have,

[TABLE]

By dominated convergence theorem we get the continuity of $g_{y}(t)$. Observe that if instead of considering the process $\eta_{t}$ at times $t=0,1,2,\dots$ we would have been considered it at times $t=0,\delta,2\delta,\dots$ we would have been obtained exactly the same result. That is, $g_{y}(n\delta)\to 0$ as $n\to\infty$, for every $\delta>0$. By [Kin63a, Theorem 2], we have the convergence in $t$,

[TABLE]

proving the result (1.3). ∎

Proof of Corollary 3.6.

We are going to prove that $\mu$ is bounded. So first, consider the process with rates $\tilde{A}$ as in the proof of Proposition 3.4 and take $V(x)=x$. For this process we have computed in (3.7)

[TABLE]

by (b). By (a) and (c) we have,

[TABLE]

Hence, $\tilde{A}$ and $A$ are $R$-positive and Yaglom limit exists. Next, observe that (a), (b) and (c) implies that the process with rates $Q$ “comes down from infinity”, meaning that $\sup_{x}\mathbb{E}_{x}(\tau_{1})<\infty$, [BMR16, Lemma 2.2 and Proposition 2.2].

We need to find a bound for $\mu$. Since $A$ is $R$-positive, $\mu$ has the following characterization

[TABLE]

To bound $\mathbb{E}_{x}|\eta_{t}|$ observe that if we start with one particle at $x$ (that we call the ancestral particle), this particles will produce during its whole life (before being absorbed) a random number of particles that can be bounded by the number of occurrences in a one-dimensional Poisson process with rate $\sup_{x}\beta(x)$ at an independent random time $T_{x}$ (the absorption time of the process with rates $Q$, started at $x$). We have just proved that the expectation of $T_{x}$ is bounded above by some constant $\kappa$. Each particle produced by the ancestral particle is born at $1$ and hence, its expected number of descendants by time $t$ is bounded by $\mathbb{E}_{1}|\eta_{t}|$. Plugging into (3.9) we get

[TABLE]

The rest of the proof is as in case 2 since we have shown $R-$positivity of $\tilde{A}$, existence of Yalgom limit for $\tilde{A}$ and boundedness of $\mu$. ∎

Acknowledgments

We thank Pablo Ferrari for several years of enlightening conversations on this topic. P. Groisman and K. Ravishankar were supported by Simons collaboration grant number 281207 awarded to K. Ravishankar. P. Groisman and A. Ferrari are partially founded by UBACYT 20020160100147BA and PICT 2015-3154. The authors want to thank NYU-Abu Dhabi where part of this work was done during P. Groisman’s visit for hospitality and support.

∎