Modelling Silicosis: The Structure of Equilibria

TL;DR

This paper analyzes the equilibrium structures of a complex coagulation-fragmentation-death model for silicosis, providing exact results, existence conditions, and asymptotic analysis for different coefficient cases.

Contribution

It offers new mathematical insights into the equilibrium multiplicity, existence, and asymptotics of the silicosis model with various coefficient configurations.

Findings

Exact multiplicity results for piecewise-constant coefficients

Conditions for existence and non-existence of equilibria

Asymptotic behavior for power law coefficient cases

Abstract

We analyse the structure of equilibria of a coagulation-fragmentation-death model of silicosis. We present exact multiplicity results in the particular case of piecewise-constant coefficients, results on existence and non-existence of equilibria in the general case, as well as precise asymptotics for the infinite series that arise in the case of power law coefficients.

Equations

$$\begin{array}[]{l}\dfrac{dM_{0}}{dt}~{}=~{}r-k_{0}xM_{0}-(p_{0}+q_{0})M_{0},\\ \mbox{}\\ \dfrac{dM_{i}}{dt}~{}=~{}k_{i-1}xM_{i-1}-k_{i}xM_{i}-(p_{i}+q_{i})M_{i},\;\;i\geq 1,\end{array}$$

$$\frac{dx}{dt}=\alpha -x\sum _{i=0}^{\infty }{k}_i{M}_i+\sum _{i=0}^{\infty }{q}_{i}i{M}_i.$$

$$M_0=\frac{r}{k_{0}x+p_0+q_0}.$$

$$M_1=\frac{r{k}_{0}x}{(k_{0}x+p_0+q_0)(k_{1}x+p_1+q_1)}.$$

$$M_i=\frac{r{x}^i{\prod }_{j=0}^{i-1}{k}_j}{{\prod }_{j=0}^i(k_{j}x+p_j+q_j)}.$$

$$M_i=\frac{r{x}^i}{k_i{\prod }_{j=0}^i(x+d_j)},i\ge 0.$$

$$p_i=\{\begin{array}{ll}1& \text{if i≤N,}\\ 0& \text{if i≥N+1,}\end{array}\text{and}{q}_i=\{\begin{array}{ll}0& \text{if i≤N,}\\ 1& \text{if i≥N+1.}\end{array}$$

$$\sum _{i=0}^{\infty }{k}_i{M}_i=k\sum _{i=0}^{\infty }{M}_i=\frac{r}{x+1/k}\sum _{i=0}^{\infty }{\left(\frac{x}{x+\text{/}k}\right)}^i=rk,$$

$$\sum _{i=0}^{\infty }i{q}_i{M}_i$$

$$\frac{\alpha }{r}-{\mathcal{F}}_{N,k}(x)=0$$

$${\mathcal{F}}_{N,k}(x):=kx\left(1-{\left(\frac{x}{x+1/k}\right)}^{N+1}\right)-(N+1){\left(\frac{x}{x+1/k}\right)}^{N+1}.$$

$$\underset{x\to +\infty }{\lim }{\mathcal{F}}_{N,k}(x)$$

$${\tilde{\mathcal{F}}}_{N,k}(y):={\mathcal{F}}_{N,k}(x(y))$$

$${\tilde{\mathcal{F}}}_{N,k}^{\prime }(y)$$

$$\sum _{i=0}^{\infty }{k}_i{M}_i<\infty \text{and}\sum _{i=0}^{\infty }i{q}_i{M}_i<\infty .$$

$$\frac{k_{i+1}{M}_{i+1}}{k_i{M}_i}=\frac{x}{x+d_{i+1}}\le \frac{x}{x+z}<1,$$

$$\frac{(i+1)q_{i+1}{M}_{i+1}}{i{q}_i{M}_i}=\frac{1+i}{i}\frac{(q_{i+1}/k_{i+1})}{q_i/k_i}\frac{x}{x+d_{i+1}}.$$

$$\frac{1+i}{i}\frac{(q_{i+1}/k_{i+1})}{q_i/k_i}\le (1+ϵ)$$

$$\frac{(i+1)q_{i+1}{M}_{i+1}}{i{q}_i{M}_i}\le \frac{1}{2}\left(\frac{x}{x+z}+1\right)<1.$$

$$\frac{\alpha }{r}=\frac{x}{x+d_0}+\frac{1}{x+d_0}\sum _{i=1}^{\infty }\left(x-i\frac{q_i}{k_i}\right)\prod _{j=1}^i\frac{x}{x+d_j}=:\mathcal{F}(x).$$

$$\mathcal{F}(x)$$

$$\frac{1}{x}\sum_{i=1}^{\infty}\Bigl{(}id_{i}-x\Bigr{)}\prod_{j=1}^{i}\frac{x}{x+d_{j}}$$

$$a_n=1-S_n.$$

$$a_n=\frac{(n+1)x^n}{{\prod }_{j=1}^n(x+d_j)}.$$

$$a_1=1-S_1=1-\frac{1}{x}\sum _{i=1}^1\left(i{d}_i-x\right)\prod _{j=1}^i\frac{x}{x+d_j}1-\frac{1}{x}(d_1-x)\frac{x}{x+d_1}=\frac{2x}{x+d_1},$$

$$a_{n+1}$$

$$0<a_{n}=(n+1)\frac{x^{n}}{\displaystyle{\prod_{j=1}^{n}(x+d_{j})}}<(n+1)\Bigl{(}\frac{x}{x+z}\Bigr{)}^{n}\longrightarrow 0,\hbox{ as }n\to+\infty,$$

$$\mathcal{F}(x)=\frac{1}{x+d_0}\sum _{i=1}^{\infty }i{\rho }_i\prod _{j=1}^i\frac{x}{x+d_j}=:\frac{1}{x+d_0}\mathcal{H}(x),$$

$$\sum _{i=1}^{\infty }{d}_i\prod _{j=1}^i\frac{x}{x+d_j}=x.$$

$$b_n=x-\sum _{i=1}^n{d}_i\prod _{j=1}^i\frac{x}{x+d_j}.$$

Full text

Research partially supported by FCT/Portugal through project UID/MAT/04459/2013.

Research partially supported by the Austrian Science Foundation FWF, project F 50-02.

1. Introduction

We examine the equations considered in [4] for the dynamics of alveolar macrophages faced with an inhalation of quartz particles in the lungs. The equations are of coagulation type, though not understood in the standard way. Thus the interest of the problem is triple: it is of medical and environmental health importance to understand how the system reacts to continuous exposure to quartz; the model shows how versatile coagulation equations are; and the mathematical structure of the system (which we do not a priori truncate as is done in [4]) presents interesting challenges. In this paper we only discuss the model itself and the structure of its equilibria, leaving proofs of global existence and stabilisation results to future work.

If we denote by $M_{i}$ the concentration of macrophages containing $i$ quartz particles (which we will call the $i$-th cohort), by $x$ the concentration of quartz, by $r$ (which can be a function of $x$ the rate of supply of new macrophages, following [4] we obtain the following equations:

[TABLE]

where $k_{i}$ is the rate of phagocytosis of a macrophage containing $i$ particles of quartz, $p_{i}$, is the transfer rate of macrophages in the $i$-th cohort to the muco-ciliary escalator, i.e. the rate of their removal together with their quartz baggage, and $q_{i}$ is the rate of death of the $i$-th cohort which results in the release of the quartz burden.

Note that unlike [4] we do not impose an upper limit on the number of particles a macrophage can contain. What is not done in [4] is to provide an equation for the evolution of the concentration of $x$; their interest is in the system dynamics following an instance of inhalation, while we are more concerned with analysing system behaviour under continuous influx of quartz. Thus we add to (1) the following equation:

[TABLE]

Here $\alpha$ the rate of inhalation of quartz.

Thus the object of our study is the system (1)–(2), considered as an infinite-dimensional dynamical system on a suitable sequence space. Before we analyse it, let us remark that hence (1)–(2) is an example of a coagulation-death system, in which the “monomers” (quartz particles) are structurally different from “clusters” (cells containing these particles); this shows the versatility of coagulation-fragmentation framework, and in particular its suitability to describe phagocytosis phenomena (e.g. of neutrophils consuming bacteria).

As in [4] we make the assumptions that $k_{i}$ and $p_{i}$ are non-increasing in $i$. We allow $q_{i}$ to grow with $i$.

The model of [4] is biologically sophisticated, also involving neutrophils and communication between neutrophils and macrophages. In (1), $r$ should express the amount of “distress” in the system, embodied in the number of macrophages with more than a sublethal load of quartz, i.e. those that are more likely to die and release their load than to be removed via the muco-ciliatory escalator. In other words, if the sublethal load is $s$ particles per cell, a biologically reasonable assumption is that $r$ is a bounded increasing function of $\sum_{i=s+1}^{\infty}M_{i}$ (see eqs. 7–8 in [4]). In the present work we take $r$ to be a constant, but our analysis here illuminates the more general case described above as well.

A simple instance of allowable coefficients for which the structure of equilibria can be analysed explicitly will be considered in section 2.1 below. The structure of the equilibria in a more general case, where the coefficients satisfy some power law relations, will be considered in section 2.2.

2. Equilibria

We start by proceeding formally and then justify our steps in the sections below. Suppose system (1)-(2) has an equilibrium. Then the $M_{0}$ equation at equilibrium can be solved for $M_{0}$ in terms of $x$ (and $r$) to give

[TABLE]

Similarly, $M_{1}$ will be given by

[TABLE]

Continuing recursively, we have

[TABLE]

Setting $d_{i}=(p_{j}+q_{j})/k_{j}$, we can rewrite this as

[TABLE]

2.1. A piecewise-constant class of coefficients

A simple instance of allowable coefficients is to take all $k_{i}$ equal to $k$ and,

[TABLE]

Then $d_{j}=1/k$, and using (3) we easily compute

[TABLE]

and

[TABLE]

Thus, plugging (4) and (5) into equation for the equilibrium quartz concentration we obtain

[TABLE]

where

[TABLE]

Proposition 1**.**

For all $r,k>0$ and $N\in{\mathbb{N}},$ there is $\alpha^{*}$ such that (6) has no solutions if $\alpha>\alpha^{*}.$

Proof. It suffices to observe that $\mathcal{F}_{N,k}(0)=0$, $\mathcal{F}^{\prime}_{N,k}(0)=k>0$ and

[TABLE]

This implies that $\mathcal{F}_{N,k}$ has an absolute maximum in ${\mathbb{R}}^{+}.$ Defining $\alpha^{*}:=\displaystyle{r\max_{{\mathbb{R}}^{+}}\mathcal{F}_{N,k}},$ the result follows. $\square$

We now prove that, for each $\alpha\in(0,\alpha^{*})$, there are exactly two solutions of (6).

Proposition 2**.**

Let $r,k>0$ and $N\in{\mathbb{N}}.$ Let $\alpha^{*}:=\displaystyle{r\max_{{\mathbb{R}}^{+}}\mathcal{F}_{N,k}}.$ Then, for every $\alpha\in(0,\alpha^{*})$ there are exactly two solutions of (6).

Proof. To prove the result we establish that $\mathcal{F}_{N,k}$ has a single stationary point in ${\mathbb{R}}^{+}$, which, then, must be the absolute maximum whose existence was established above. This, together with the already proved facts that $\mathcal{F}_{N,k}(0)=0$ and $\lim_{x\to+\infty}\mathcal{F}_{N,k}(x)=0$, proves the result.

Let $y:=\frac{x}{x+1/k}$. Then $x=\frac{y/k}{1-y}$ and

[TABLE]

Since $\frac{dx}{dy}=\frac{1}{k}\frac{1}{(1-y)^{2}}>0$, we have $\operatorname{sgn}\widetilde{\mathcal{F}}^{\prime}_{N,k}(y)=\operatorname{sgn}\mathcal{F}^{\prime}_{N,k}(x(y))$. Thus we need only to study the function in the new variable $y\in[0,1)$. Observing that $f^{\prime}_{N}(y)=-N(N+1)y^{N-1}+N(N+1)y^{N}=-N(N+1)y^{N-1}(1-y),$ we have

[TABLE]

Let us consider the polynomial $p_{N}$ in $[0,1].$ It is clear that $p_{N}(0)=1$ and $p_{N}(1)=0.$ Its derivative is $p_{N}^{\prime}(y)=N(N+1)y^{N-1}q_{N}(y)$, where $q_{N}(y):=-(N+1)+(2N+3)y-(N+2)y^{2}.$ We easily conclude that the zeros of $q_{N}(y)$ are $y_{1}=\frac{N+1}{N+2}$ and $y_{2}=1,$ and that $\operatorname{sgn}\left(y-\frac{N+1}{N+2}\right)\operatorname{sgn}q_{N}(y)>0.$ This means that $p_{N}$ has a minimum at $y=\frac{N+1}{N+2}$ and must be an increasing function in the interval $\left(\frac{N+1}{N+2},1\right)$. Since $p_{N}(1)=0,$ this implies the value of $p_{N}(y)$ at $y=\frac{N+1}{N+2}$ must be negative, which, together with $p_{N}(0)>0$ and the fact that $p_{N}$ is strictly decreasing in $\left(0,\frac{N+1}{N+2}\right)$, means, due to the intermediate value theorem, that there is one, and only one, zero of $p_{N}$ in this set, and hence in $(0,1)$, i.e., there is a single stationary point of $\mathcal{F}_{N,k}$ in ${\mathbb{R}}^{+}$. $\square$

2.2. Power type coefficients

We consider now the more complex case of coefficients satisfying some power relations.

Theorem 3**.**

Let $M_{i}$ be given by (3). Assume that $z=\inf_{i}d_{i}>0$. Assume also that $q_{i}/k_{i}$ grows no faster than a power of $i$. Then for all $x\geq 0$

[TABLE]

Proof. This follows by the Ratio Test, as

[TABLE]

for all $i$. Also,

[TABLE]

Pick $\epsilon={\displaystyle\frac{1}{2}\left(\frac{z+x}{x}-1\right)}$. We can find $N=N(x)$ such that

[TABLE]

for all $i\geq N$. But then for all $i\geq N$ we have that

[TABLE]

$\square$

So the equation for equilibrium quartz concentration can be written in the form

[TABLE]

Our first main result provides a quite general sufficient condition for the existence of equilibria.

Theorem 4**.**

Assume that $z:=\inf_{i}d_{i}>0$ and also that $q_{i}/k_{i}$ grows no faster than a power of $i$. Let $\rho_{i}:=p_{i}/k_{i}.$

If $d_{i}=o(i\rho_{i})$ (as $i\to\infty$) then $\mathcal{F}(x)\to\infty$ (as $x\to\infty$) and consequently we have an equilibrium (8) for all $\alpha,r$.

If $i\rho_{i}=O(d_{i})$ (as $i\to\infty$) then $\mathcal{F}(x)$ is bounded. Thus there exists $m>0$ such that there exists an equilibrium for $\alpha/r<m$ and no equilibrium for $\alpha/r\geq m$.

Finally if $i\rho_{i}=o(d_{i})$ (as $i\to\infty$) then $\mathcal{F}(x)$ is bounded and we have $\mathcal{F}(x)\to 0$ (as $x\to\infty$). In this case there exists $m>0$ such that there exists an equilibrium for $\alpha/r\leq m$ and no equilibrium for $\alpha/r>m$.

In order to prove Theorem 4 we have to study $\mathcal{F}(x)$ in more detail.

Since $d_{i}:=\frac{p_{i}+q_{i}}{k_{i}}=\rho_{i}+\frac{q_{i}}{k_{i}}$, we have

[TABLE]

Proposition 5**.**

With the above assumptions and notation, we have that $\mathcal{G}(x)=0,$ $\forall x\geq 0$.

Proof. Since $\mathcal{G}(0)=0$ it is sufficient to consider the case $x>0$. We prove that the sum $S$ of the series

[TABLE]

is equal to $1$ for all values of $x>0.$ Let $S_{n}$ denote the partial sums of (9) and set

[TABLE]

We will show by induction that

[TABLE]

Obviously we have for $n=1$

[TABLE]

as required. Assume (10) is true for some $n$. Then

[TABLE]

This proves (10) for all $n\geq 1$.

Now $a_{n}\to 0$ as $n\to\infty$ follows trivially from

[TABLE]

So we conclude that $S=\lim S_{n}=\lim(1-a_{n})=1$. Hence $\mathcal{G}(x)=0$ holds also for all $x>0.$ $\square$

Using Proposition 5 we conclude that, with the power law assumptions on the coefficients, $\mathcal{F}$ can be written as

[TABLE]

With the help of the next proposition we can get some information on the growth order of $\mathcal{H}(x)$.

Proposition 6**.**

Assume that $z=\inf_{i}d_{i}>0$. Then we have, for $x\geq 0$,

[TABLE]

Proof. The equality is trivially satisfied for $x=0$. Thus we just have to consider the case $x>0$, where we set

[TABLE]

We prove by induction that

[TABLE]

This is clearly true for $n=1$:

[TABLE]

Now assume that (12) is satisfied for some $n\geq 1$. Then we have

[TABLE]

as proposed. Finally we have

[TABLE]

This implies $\lim b_{n}=0$ and proves the proposition for $x>0$. $\square$

We are now ready to prove Theorem 4.

Proof. First suppose that $d_{i}=o(i\rho_{i})$ (as $i\to\infty$). Fix some $\varepsilon>0$ and suppose that $d_{i}\leq\varepsilon i\rho_{i}$ for $i\geq i_{0}=i_{0}(\varepsilon)$. Then we have (also by applying Proposition 6)

[TABLE]

Consequently

[TABLE]

Since $\varepsilon>0$ can be arbitrarily chosen we have $\mathcal{H}(x)/x\to\infty$ (as $x\to\infty$) and, thus,

[TABLE]

Since $\mathcal{F}(0)=0$ and $\mathcal{F}(x)$ is continuous it follows that there exists an equilibrium in all cases.

Next suppose that $i\rho_{i}=O(d_{i})$, that is, there is a constant $K>0$ such that $i\rho_{i}\leq Kd_{i}$ for all $i\geq 1$. Hence,

[TABLE]

and consequently $\mathcal{F}(x)$ is bounded. Clearly, if we set

[TABLE]

then there exists no equilibrium if $\alpha/r>m$ and again since $\mathcal{F}(0)=0$ and by continuity there is an equilibrium if $\alpha/r<m$.

Finally if $i\rho_{i}=o(d_{i})$ (as $i\to\infty$) then is follows (as above) that $F(x)=o(x)$ (as $x\to\infty$). By continuity there exists

[TABLE]

Hence, then there exists no equilibrium if $\alpha/r>m$ and an equilibrium if $\alpha/r\leq m$. $\square$

Consider now the case where the coefficients are given by the following power laws:

[TABLE]

for $i\in{\mathbb{N}}^{+}$ and nonnegative constants $p,q$ and $k$, Let $p_{0},q_{0}$ and $k_{0}\neq 0$ be given. Then, writing $a:=q+k\geq 0$ and $b:=k-p\in{\mathbb{R}},$ we have $d_{i}=i^{a}+i^{b}$ and $\rho_{i}=i^{b}.$

By a direct application of Theorem 4 we get the following property:

Corollary 7**.**

Suppose that $a\geq 0$ and $b>-2$.

If $b>a-1$, (1)–(2) has an equilibrium for all $\alpha/r$.

If $b=a-1$, there is a value $m>0$ such that for $\alpha/r<m$, (1)–(2) has an equilibrium and no equilibria if $\alpha/r\geq m$.

If $b<a-1$, there is a value $m>0$ such that for $\alpha/r\leq m$, (1)–(2) has an equilibrium and if $\alpha/r>m$, there are none.

3. Precise Asymptotics

It is also an interesting problem to obtain precise asymptotics for the case where $d_{i}=i^{a}+i^{b}$ and $\rho_{i}=i^{b}$, $b<a$. In order to make our analysis slightly easier we will concentrate on the case

[TABLE]

We will use the following notation:

[TABLE]

Then we have the following result:

Theorem 8**.**

Suppose that $a>0$, $b>-2$ and let $\mathcal{K}_{a,b}(x)$ be given by (13). Then as $x\to\infty$, $\mathcal{K}_{a,b}(x)$ admits an asymptotic expansion such that

[TABLE]

if $b>-1$ and

[TABLE]

if $b\leq-1$.

The proof is mainly based on the following asymptotic series representation:

Lemma 9**.**

Let $R(a,A;v)$ denote the infinite sums

[TABLE]

where $a,v>0$, and $A$ is real. The following holds:

If $(A+1)/(a+1)$ is different from $0,-1,-2,\ldots$, then as $v\to 0$,

[TABLE]

where $\zeta(\cdot)$ is the Riemann zeta function.

• 2.

If $(A+1)/(a+1)=-k_{0}$ for some integer $k_{0}\geq 0$ then, as $v\to 0$,

[TABLE]

where $H_{k}=1+\frac{1}{2}+\cdots\frac{1}{k}$ denotes the $k$-th harmonic number, $H_{0}=0$, and $\gamma=0.5772156\ldots$ is the Euler–Mascheroni constant.

Proof. We recall (see, e.g., [1, Part I]) that the Mellin transform of a function $f(v)$ is given by

[TABLE]

and converges usually in a strip $a_{1}<\Re(s)<a_{2}$. Under suitable regularity assumptions (for example that $f(v)$ is continuous and of bounded variation) the function $f(v)$ can be recovered from the integral

[TABLE]

where $a_{1}<C<a_{2}$.

In our case it is an easy exercise to show that the Mellin transform of $R(a,A;v)$ is given by

[TABLE]

The integral converges for \Re(s)>\max\bigl{(}0,\frac{A+1}{a+1}\bigr{)}. Consequently we have

[TABLE]

where C>\max\bigl{(}0,\frac{A+1}{a+1}\bigr{)}. Since the $\Gamma$-function decreases exponentially fast on vertical lines we could replace the limit $\lim_{T\to\infty}\int_{C-iT}^{C+iT}$ by the indefinite integral $\int_{C-i\infty}^{C+i\infty}$.

The function $\Gamma(s)\zeta((a+1)s-A)$ has a meromorphic continuation to the whole complex plane. The only singularities are simple poles coming from $\Gamma(s)$ at $s=-k$ with residue

[TABLE]

and the simple pole of $\zeta((a+1)s-A)$ at $s=\frac{A+1}{a+1}$ with residue

[TABLE]

The idea is to shift the integral in (16) to the left and to collect residues of the polar singularities that are passed. There are (again) no convergence problems of the integral due to the $\Gamma$ factor.

Assume first that $\frac{A+1}{a+1}$ is different from $0,-1,-2,\ldots$. If we shift the integral to $\Re(s)=-M-\frac{1}{2}$, where $M$ is a positive integer and $-M<\frac{A+1}{a+1}$ then we have

[TABLE]

which proves the first part of the lemma.

If $(A+1)/(a+1)=-k_{0}$ for some integer $k_{0}\geq 0$ then $\Gamma(s)$ and $\zeta((a+1)s-A)$ create a double pole at $s=-k_{0}$ with residue

[TABLE]

of the resulting function. This explains the difference from the first case and completes the proof of the lemma. $\square$

We also need representations for finite sums of powers of integers that can be deduced from the Euler–McLaurin summation formula, see e.g. [2, Chapter 9].

Lemma 10**.**

We have the following representations or asymptotic series representation, resp., for the sums $\sum_{j=1}^{n}j^{a}$, $a>0$:

1. a)

If $a$ is a non-negative integer,

[TABLE]

where $B_{2k}$ are the Bernoulli numbers. 2. b)

If $a$ is a real number different from the non-negative integers, we have the asymptotic series expansion

[TABLE]

We are now ready to prove Theorem 8. Let us write $\mathcal{K}_{a,b}(x)=\sum_{i=1}^{\infty}P_{i}(x)$, where $P_{i}(x):=i^{b+1}\prod_{j=1}^{i}\frac{1}{1+j^{a}/x}$, and note that, as $x\to\infty$,

[TABLE]

First of all we prove that $P_{i}(x)$ do not contribute significantly to $\mathcal{K}_{a,b}(x)$ if $i>x^{3/(3a+2)}$.

Lemma 11**.**

We have

[TABLE]

for some constants $c,K>0$.

Proof. First, we assume that $i^{a}>\alpha x$, where $\alpha$ will be chosen later. Then we have

[TABLE]

Since $i!\geq(i/e)^{i}$ and $i^{a}>\alpha x$ we thus obtain

[TABLE]

Consequently, if we choose $\alpha=e^{2a+(1+b)}$, we have that $P_{i}(x)\leq e^{-ai}$ and hence

[TABLE]

for some constants $c_{2},K_{1}>0$.

We now assume that $x^{3/(3a+2)}<i\leq\alpha^{1/a}x^{1/a}$, with $\alpha$ chosen as above. In this case we have that $j^{a}/x\leq\alpha$ so that there exists a constant $c_{3}>0$ such that

[TABLE]

for all $j\leq i$. Consequently there exists a constant $c_{4}>0$ such that

[TABLE]

Note that for every real $\rho$, every $\sigma>0$ and $\kappa_{1},\kappa_{2}$ such that $0<\kappa_{1}<\kappa_{2}$ there is a constant $D$ depending on these four numbers such that for all $y>0$

[TABLE]

Hence there are constants $K_{2}>0$ and $c_{5},c_{6}>0$ such that

[TABLE]

Now pick $c=\min\{c_{2},c_{6}\}$ and $K=\max\{K_{1},K_{2}\}$ to complete the proof of the lemma. $\square$

Thus, it remains to consider $i$ with $i\leq x^{3/(3a+2)}$. In this case we certainly have $j^{a}/x\to 0$ as $x\rightarrow\infty$, so we can use the Taylor expansion of $\log(1+z)$ to proceed further. From this Taylor series expansion, it follows that for every $L\geq 1$ we have uniformly for $j\leq i$

[TABLE]

and consequently

[TABLE]

In order to handle these terms we will use Lemma 10.

With the help of the representation (19) and Lemma 10 we see that $\sum_{j=1}^{i}\log\left(1+j^{a\ell}/x\right)$ is dominated by $i^{a\ell+1}/((a\ell+1)x)$ followed by smaller order terms. Here we have to distinguish between $\ell=1$ and $\ell\geq 2$. For $\ell=1$ the dominating term $i^{a+1}/((a+1)x)$ is unbounded if $i\leq x^{3/(3a+2)}$ whereas the next order term $i^{a}/x$ (and all following terms) are bounded (in order) by $x^{-2/(3a+2)}$. So all of them go to zero if $x\to\infty$. It $\ell\geq 2$ then the dominant terms $i^{a\ell+1}/((a\ell+1)x)$ (and, thus, all following terms) will go to zero, too. They are bounded (in order) by $x^{(3-2\ell)/(3a+2)}\leq x^{-1/(3a+1)}$.

Summing up, we obtain for $i\leq x^{3/(3a+2)}$

[TABLE]

where $\tilde{S}$ collects terms of the form $\log\frac{1}{x}$ that go to zero. Hence by using the Taylor series of the exponential function we have

[TABLE]

which leads again to an asymptotic series representation for $P_{i}(x)$ of the form

[TABLE]

where $\tilde{T}$ collects terms of the form $\hbox{const}\cdot i^{A}/x^{B}$ (with real $A$ and integer $B\geq 1$) that go to zero if $i\leq x^{3/(3a+2)}$.

This discussion shows that we are finally led to consider sums of the form

[TABLE]

Since the sum of the missing terms can be estimated by

[TABLE]

for some constant $c_{7}>0$, it is sufficient to consider infinite sums of the form analysed in Lemma 9.

We are now ready to complete the proof of Theorem 8.

Since (20) are asymptotic series for $P_{i}(x)$, it follows that we can consider them always as finite sums plus an error term of the same form. Thus we can sum over them, at least for $i\leq x^{3/(3a+2)}$. However, by Lemma 11 we can extend this summation over all $i\geq 1$ since the resulting error is negligible.

Considering the terms in $\tilde{T}$, we see that the asymptotic representation is different for the case of $a$ being a positive integer, and for non-integer real positive $a$.

In the first case, by applying the (17) of Lemma 10 and observing that we just get positive powers of $n$ in the representation of the sums $\sum_{j=1}^{i}j^{\ell a}$, $1\leq\ell<L$, the asymptotic series expansion (20) of $P_{i}(x)$, $i\leq x^{3/(3a+2)}$ can be written in the form

[TABLE]

where $A,B$ are now integers and $\tilde{c}_{A,B}$ are real constants.

This means that we also get an asymptotic series representation of $\mathcal{K}_{a,b}(x)$ of the form

[TABLE]

For non-integer $a$ we we can proceed in the same way as in the integer case. There are, however, some differences in the course of the computations. First of all the sums $\sum_{j=1}^{i}j^{\ell a}$ do not have an explicit representations. By (18) of Lemma 10, we obtain an asymptotic series expansion that contains also negative powers of $i$, namely $i^{a+1-k}$ for any $k\geq 0$. This leads to an asymptotic series expansion for $P_{i}(x)$ of the form

[TABLE]

where $c_{A,B}$ are again real constants and the sum ranges over all (even negative) integers $A\leq\lfloor B/2\rfloor$. In completely the same way as above we get from that an asymptotic series expansion for $\mathcal{K}_{a,b}(x)$:

[TABLE]

Now we are going to use the information in Lemma 9 to understand the leading terms and the order of the remainder in (21) and (22).

In both cases, of integer and non-integer $a$, we have that for any fixed integer $B$, the expression ${\displaystyle\frac{b+2+A}{a+1}-B}$ (for integer $a$) and ${\displaystyle\frac{b+2+aB+A}{a+1}-B}$ (for non-integer $a$) are maximised by taking the largest allowable $A$ in these two cases to give

[TABLE]

and this is maximised by picking $B=1,2$ to give $(b+1)/(a+1)$. Hence if $b>-1$, the asymptotic series (21) and (22) give, using (14)

[TABLE]

If $b\leq-1$, we get

[TABLE]

with the contributions to the $O(1)$ term coming from $R(a,1+b;1/((a+1)x))$ term if $b<-1$. If $b=-1$, the contributions to the $O(1)$ term come again from the $R(a,1+b;1/((a+1)x))$ term, and, if $a$ is an integer, from the $B=1,A=1$ term in (21); if $a$ is not an integer, the additional contributions come from the $B=1,A=0$ and the $B=2,A=1$ terms in (22).

Remarks: 1. Much more can be said using Lemma 9. For example, if $b=-2$ the leading term of the asymptotic expansion of $K_{a,b}(x)$ will be of order $\log x$.

2. The case $b=-1$, $a=1$ can be solved explicitly. There we have

[TABLE]

3. It is easy to show that

[TABLE]

Another case that can be solved explicitly is the case of $\mathcal{H}_{1,0}(x)$,

[TABLE]

In this case MAPLE can compute the series and its asymptotics. We have

[TABLE]

where

[TABLE]

is the incomplete $\Gamma$-function. The asymptotic expansion as $x\rightarrow\infty$ is

[TABLE]

4. Further Remarks and Conclusions

Corollary 7 is clearly non-optimal. Numerical evidence shows that if $b\geq a-1$ the equilibrium is unique for all $\alpha/r$ and there are at most two equilibria if $b<a-1$. To prove such a result one might try to use the machinery of Pinelis [3]. The result for $b\geq a-1$ would follow if one could prove that $\mathcal{H}(x)$ is convex and for $b<a-1$, by showing $\mathcal{H}(x)$ is concave. However, these results are difficult to obtain even for $\mathcal{K}_{a,b}(x)$. It does seem that for each fixed $x$, the second derivative of $\mathcal{K}_{a,b}(x)$ is an increasing function of $b$.

In the case of $\mathcal{H}_{1,0}(x)$ considered above, numerically $\mathcal{H}_{1,0}$ appears to be convex. Hence if we could prove that, and monotonicity of the second derivative in $b$ for $\mathcal{H}_{a,b}(x)$ for every fixed $x$, the desired results for $\mathcal{H}_{a,b}(x)$ would follow by the argument of Pinelis [3].

If $b<a-1$, a possible strategy for proving, for example, that $\mathcal{K}_{a,b}(x)$ is concave is to consider partial sums of the infinite sum. For $b\geq a-1$ this strategy also works but is more interesting because $x/(x+1)$ is concave. In that case, if $S_{n}(x)$ is the partial sum, numerics indicate that $S_{n}^{\prime\prime}(x)$ is positive for $0<x<x_{n}$ and $x_{n}\rightarrow\infty$ s $n\rightarrow\infty$.

Note that for all integer values of $a,b$ MAPLE can compute $\mathcal{F}_{a,b}(x)$ in terms of hypergeometric functions. This, however, does not seem very useful.

From the biological point of view, $b$ measures the efficiency of the muco-ciliatory escalator, while $a$ measures its inefficiency due to release of quartz in the lungs by macrophages with supercritical load. Our results show that the ratio $(b+2)/(a+1)$ is crucial in establishing whether the system can deal with the quartz load; if it is less or equal to 1, there is a deposition rate $\alpha$ that will overwhelm it, no matter what $r$ is.

In summary, we have completed the model of [4] by including an equation for the evolution of quartz concentration. The resulting mathematical object is a challenging system of coagulation–death equations that requires non-trivial asymptotic ideas in the discussion of the structure of equilibria. Of course the analysis in the paper is only part of the necessary mathematical work; one also needs to establish global existence (using finite-dimensional truncations or methods of semigroup theory) and stabilisation to equilibria (for example, by exhibiting a suitable Lyapunov function).

Acknowledgements

FPdC and MG would like to thank P. Freitas and D. Pritchard for valuable discussions. FPdC acknowledges financial support provided by the University of Strathclyde David Anderson Research Professorship, and all the authors are grateful to D. Zeilberger for helping to set up the collaboration.