# Global existence and asymptotic behavior for a reaction-diffusion system
with unbounded coefficients

Mohamed Majdoub and Nasser-eddine Tatar  Department of Mathematics, College of
Science, Imam Abdulrahman Bin Faisal University, P. O. Box 1982, Dammam, Saudi
Arabia.
Basic and Applied Scientific Research Center, Imam Abdulrahman Bin Faisal
University, P.O. Box 1982, 31441, Dammam, Saudi Arabia<EMAIL_ADDRESS>King Fahd University of Petroleum and Minerals, Department of Mathematics, IRC
for Intelligent Manufacturing and Robotics, Dhahran, 31261, Saudi Arabia.
<EMAIL_ADDRESS>
###### Abstract.

We consider a reaction-diffusion system which may serve as a model for a
ferment catalytic reaction in chemistry. The model consists of a system of
reaction diffusion equations with unbounded time dependent coefficients and
different polynomial reaction terms. An exponential decay of the globally
bounded solutions is proved. The key tool of the proofs are properties of
analytic semigroups and some inequalities.

###### Key words and phrases:

Analytic semi-group, exponential decay, fractional operator, reaction-
diffusion system, sectorial operator.

###### 2010 Mathematics Subject Classification:

35K57, 35A01, 35K57,35B40, 35Q92

##  1\. Introduction

In the past three decades a strong interest in the questions of global
existence and asymptotic behavior of solutions to various classes of reaction-
diffusion systems has been witnessed. A particular attention was paid to
autonomous systems and different sufficient conditions were assumed to ensure
global existence of classical solutions. A main difficulty in proving global
existence is the lack of maximum principle estimates or invariant regions.
See, among many, [1, 4, 3, 10, 12, 14, 15, 17]. Recently, the global existence
of classical solutions was shown in [2] for quasi-positive nonlinearities
having a (slightly super-) quadratic growth and obeying a mass control
assumption. This, in particular, includes the cases of mass conservation and
mass dissipation. See also [11] for initial data of low regularity.

Our goal in this paper is to study the following non-autonomous system

$\displaystyle\left\\{\begin{array}[]{ccl}u_{t}-(d_{1}\Delta-
b_{1})u&=&a_{1}(t)w^{m}-a_{2}(t)u^{n}v^{k},\\\ v_{t}-d_{2}\Delta
v&=&\left(a_{1}(t)+a_{3}(t)\right)w^{m}-a_{2}(t)u^{n}v^{k},\\\
w_{t}-(d_{3}\Delta-
b_{2})w&=&-\left(a_{1}(t)+a_{3}(t)\right)w^{m}-a_{4}(t)w+a_{2}(t)u^{n}v^{k},\\\
\frac{\partial u}{\partial\nu}=\frac{\partial v}{\partial\nu}=\frac{\partial
w}{\partial\nu}&=&0,\\\
(u,v,w)(x,0)&=&(u_{0}(x),v_{0}(x),w_{0}(x)),\end{array}\right.$ (1.6)

where $t>0$, $x\in\Omega\subset\mathbf{R}^{N},$ $N\geq 1$, $\Omega$ is a
bounded domain with smooth boundary $\partial\Omega$ and
$\partial/\partial\nu$ denotes the outward normal derivative on
$\partial\Omega.$ The diffusion coefficients $d_{i},$ $i=1,2,3$ are assumed
positive constants and $u_{0},$ $v_{0}$ and $w_{0}$ are given nonnegative
bounded initial functions. The coefficients $a_{i}(t)$ are of the form
$t^{\sigma_{i}}h_{i}(t),$ where $\sigma_{i}\geq 0,$ $i=1,2,3$ and
$\sigma_{4}=0.$ The functions $h_{i}(t)$ are nonnegative continuous functions.
The powers $m,$ $n$ and $k$ are assumed greater than one.

Systems with time dependent nonlinearities were considered in [7] by Kahane.
Specifically, a system of the form

$\left\\{\begin{array}[]{l}-u_{t}+Lu=f(x,t,u,v)\;in\;\Omega\times(0,\infty)\\\
-v_{t}+Mv=g(x,t,u,v)\;in\;\Omega\times(0,\infty)\end{array}\right.$

where $L$ and $M$ are uniformly elliptic operators, with boundary conditions
of Robin type is studied. He proved that the solution converges to the
stationary state, that is the solution of the limiting elliptic problem. To
this end he assumed that

$f(x,t,u,v)\rightarrow\bar{f}(x,u,v)$

and

$g(x,t,u,v)\rightarrow\bar{g}(x,u,v)$

uniformly in $\Omega$ and $(u,v)$ in any bounded subset of the first quadrant
in $\mathbf{R}^{2}.$ The matrix formed by the partial derivatives
$\bar{f}_{u},\bar{f}_{v},\bar{g}_{u}$ and $\bar{g}_{v}$ satisfies a column
diagonal dominance type condition. In the present work we do not make such
restrictions.

Let us return now to our system (1.6). In the case $m=n=k=1,$ $b_{1}=b_{2}=0$
and $a_{i}(t)\equiv a_{i}$ are constants, this problem has been studied by
Wang [19]. If moreover $a_{3}=a_{4}=0,$ then this problem may also serve as a
model for sugar transporting into red blood cells (see Ruan [18], Rothe [17],
Feng [3], Morgan [14] and references therein). In particular, Wang proved a
convergence result and an exponential decay result in $C^{\mu}(\bar{\Omega})$,
$\mu\in[0,2).$ Unfortunately, his method seems to be not valid for the present
problem (1.6) because of the unboundedness of the coefficients.

To our knowledge, a system with time-dependent coefficients like (1.6) is
still not well studied. Nevertheless, using similar ideas as in [2], we are
able to prove global existence of classical solutions to (1.6) provided that
$0\leq m\leq 2+\varepsilon$, $n,k\geq 0$ and $n+k\leq 2+\varepsilon$ for
sufficiently small $\varepsilon>0$. More precisely, we assume that

$0\leq m\leq 2+\varepsilon\quad\mbox{and}\quad n+k\leq 2+\varepsilon,$ (1.7)
$u_{0},v_{0},w_{0}\in L^{1}(\Omega)\cap L^{\infty}(\Omega)\quad\mbox{and}\quad
u_{0},v_{0},w_{0}\geq 0.$ (1.8)

######  Theorem 1.

Suppose that assumption (1.8) is fulfilled. Then, there exists $\varepsilon>0$
such that if (1.7) is satisfied, (1.6) possesses a unique classical global
solution.

Our main task here is to study the asymptotic behavior and the decay rate of
global classical solutions obtained in Theorem 1. To this end we shall adapt
the methods used in Hoshino [5] for the former task and the methods used in
Kirane and Tatar [9] for the latter goal. To state our results in a clear way,
we shall fix some notations. For $\alpha\in(0,1)$, we denote by

$q^{*}=q^{*}(\alpha)=\left\\{\begin{array}[]{l}\frac{2}{\alpha}-1,\;if\;\frac{1}{2}\leq\alpha<1\\\
2,\;if\;0<\alpha<\frac{1}{2}\end{array}\right.$

and $q$ the Lebesgue conjugate exponent of $q^{*}$. Let
$r=\frac{1-\alpha}{\alpha}$. First, we state a convergence result of the
solutions in the space $C^{\mu}(\Omega),$ $\mu\in[0,2).$ This result will be
needed in the sequel. Without loss of generality, we shall assume that
$b_{1}=b_{2}=b.$

######  Theorem 2.

Suppose that $1-q\alpha>0$, $1+q(\sigma_{i}-\alpha l)>0,$ $i=1,2,3$ for some
$l$ such that $1<l<\min\\{m,n,k\\},$ $h_{i}\in L^{q^{*}}(0,\infty),$
$i=1,2,3,4$ and $a_{1}(t)\leq Ca_{3}(t)$ for some positive constant $C$,
$\forall t>0.$ Then, for every $\mu\in[0,2)$

$u(t)\rightarrow 0,v(t)\rightarrow v_{\infty}\;and\;w(t)\rightarrow
0\;in\;C^{\mu}(\Omega)\;as\;t\rightarrow\infty$

where
$v_{\infty}=\left|\Omega\right|^{-1}\left\\{\int\limits_{\Omega}w_{0}dx-\int\limits_{0}^{\infty}\int\limits_{\Omega}(h_{4}(s)+b)w(s)dxds\right\\}.$

######  Remark 1.

The proof of Theorem 2 is similar to the proof of Theorem 2.2 in [5] with
minor modifications that will be clear from the proof of our next and main
result. It is therefore omitted.

To state our next result, we suppose that $N$ and $p$ satisfy

$2\alpha>\frac{N}{p}\quad\mbox{and}\quad\max\left\\{1,\frac{N(m-1)}{2p\alpha},\frac{N(n-1)}{2p\alpha}\right\\}<\min\\{m,n\\}.$
(1.9)

######  Theorem 3.

Assume that the hypotheses of Theorem 2 and (1.9) are fulfilled. If

$(l-1)\int\limits_{0}^{\infty}h(s)ds<\log\left(\frac{1+C_{0}^{l-1}}{C_{0}}\right)$

for some constant $l$ such that
$\max\left\\{1,\frac{N(m-1)}{2p\alpha},\frac{N(n-1)}{2p\alpha}\right\\}<l<\min\\{m,n\\}$
and a constant
$C_{0}=C_{0}(\left\|u_{0}\right\|_{p},\left\|w_{0}\right\|_{p})$ and a
function $h$ to be determined in the proof, then for $\mu\in[0,2)$ and $N,$
$p$ such that $0\leq\mu<2\alpha-\frac{N}{p}$ we have

  1. (a)

$\left\|u\right\|_{C^{\mu}(\bar{\Omega})},$
$\left\|w\right\|_{C^{\mu}(\bar{\Omega})}\leq
Ce^{-(b-\varepsilon)t}\left(\left\|u_{0}\right\|_{p}+\left\|w_{0}\right\|_{p}\right)$,
where $0<\varepsilon<b$, as $t\rightarrow\infty.$

  2. (b)

     * (i)

If $d_{2}\lambda<lb,$

$\left\|v-v_{\infty}\right\|_{C^{\mu}(\bar{\Omega})}\leq
Ce^{-\min\left\\{(b-\varepsilon),d_{2}\lambda\right\\}t}\;as\;t\rightarrow\infty.$

     * (ii)

If $d_{2}\lambda\geq lb$ and for $i=1,2,3$

$\int\limits_{0}^{t}e^{q^{*}\rho
s}h_{i}^{q^{*}}(s)ds=O(e^{q^{*}\tilde{\rho}t})\;as\;t\rightarrow\infty$

for some $\rho>d_{2}\lambda-l(b-\varepsilon)$ and $\tilde{\rho}<d_{2}\lambda,$
then

$\left\|v-v_{\infty}\right\|_{C^{\mu}(\bar{\Omega})}\leq
Ce^{-\min\left\\{(b-\varepsilon),d_{2}\lambda-\tilde{\rho}\right\\}t}\;as\;t\rightarrow\infty.$

##  2\. Preliminaries

In this section we present the notation that will be used in this paper and
prepare some material which will be useful in our proofs.

By $W^{l,p}(\Omega)$ we denote the usual Sobolev space of order $l\geq 0$ for
$1\leq p\leq\infty$. The space $C^{\mu}(\bar{\Omega}),\mu\geq 0$ is the Banach
space of $[\mu]-$times continuously differentiable functions in $\bar{\Omega}$
whose $[\mu]-$th order derivatives are Hölder continuous with exponent
$\mu-[\mu].$

######  Definition 1.

For $p\in(1,\infty)$, we define

$D(A_{p})=D(B_{p})=D(G_{p})=\left\\{y\in W^{2,p}(\Omega):\partial
y/\partial\nu\left|{}_{\partial\Omega}\right.=0\right\\},$

$A_{p}y=-(d_{1}\Delta-b_{1})y,$ $B_{p}y=-d_{2}\Delta y,$
$G_{p}y=-(d_{3}\Delta-b_{2})y.$

The operators $-A_{p},$ $-B_{p}$ and $-G_{p}$ defined in this way are
sectorial operators (see Henry [4]) and generate analytic semigroups
$\left\\{e^{-tA_{p}}\right\\}_{t\geq 0},$ $\left\\{e^{-tB_{p}}\right\\}_{t\geq
0}$ and $\left\\{e^{-tG_{p}}\right\\}_{t\geq 0},$ respectively.

######  Definition 2.

By $Q_{0}$ and $Q_{+}$ we designate the following projection operators, for
$y\in L^{p}(\Omega)$, $p\in(1,\infty)$

$Q_{0}y=\left|\Omega\right|^{-1}\int\limits_{\Omega}y(x)dx\;and\;Q_{+}y=y-Q_{0}y$

where $\left|\Omega\right|$ is the volume of $\Omega.$

######  Definition 3.

The restriction of $B_{p}$ onto $Q_{+}L^{p}(\Omega)$ will be denoted by
$B_{p+}$ i.e. $B_{p+}=B_{p}\left|{}_{Q_{+}L^{p}(\Omega)}\right..$

For $p\in(1,\infty),$ the operator $B_{p+}$ generates an analytic semigroup
denoted by $\left\\{e^{-tB_{p+}}\right\\}_{t\geq 0}$ in $Q_{+}L^{p}(\Omega).$

We define the fractional powers of the above operators in the usual way (see
Henry [4]).

######  Lemma 1.

Let $A$ be a sectorial operator in $X=L^{p}(\Omega),$ $1\leq p<\infty$ with
$D(A)=X^{1}\subset W^{m,p}(\Omega)$ for some $m\geq 1.$ Then, for
$0\leq\alpha\leq 1$, we have $X^{\alpha}\subset C^{\mu}$ when
$0\leq\mu<m\alpha-\frac{N}{p}.$

######  Lemma 2.

Let $\lambda$ denote the least positive eigenvalue of the Laplacian with
homogeneous Neumann boundary condition. Let $p\in(1,\infty)$. For every
$\alpha\in[0,1)$ there exist positive constants $C_{i},$ $i=1,2,3$ such that
for $t>0$ and $y\in L^{p}(\Omega)$

$\left\|A_{p}^{\alpha}e^{-tA_{p}}y\right\|_{p}\leq
C_{1}t^{-\alpha}e^{-b_{1}t}\left\|y\right\|_{p},$

$\left\|B_{p+}^{\alpha}e^{-tB_{p+}}Q_{+}y\right\|_{p}\leq
C_{2}t^{-\alpha}e^{-d_{2}\lambda t}\left\|Q_{+}y\right\|_{p},$

$\left\|G_{p}^{\alpha}e^{-tG_{p}}y\right\|_{p}\leq
C_{3}t^{-\alpha}e^{-b_{2}t}\left\|y\right\|_{p}.$

######  Lemma 3.

Let $p\in(1,\infty),$ $l\geq 1$ and $\alpha\in[0,1).$ Then, there exists a
positive constant $C$ such that

$\left\|y\right\|_{pl}\leq
C\left\|A_{p}^{\alpha}y\right\|_{p}^{\theta}.\left\|y\right\|_{p}^{1-\theta}$

where $\theta$ satisfies

$\frac{N(l-1)}{2pl\alpha}<\theta<1.$

See Henry [4] for the proofs of Lemmas 1-2-3.

######  Lemma 4.

Let $\alpha\in[0,1)$ and $\beta\in\mathbf{R}.$ There exists a positive
constant $C=C(\alpha,\beta)$ such that

$\int\limits_{0}^{t}s^{-\alpha}e^{\beta
s}ds\leq\left\\{\begin{array}[]{l}Ce^{\beta t}\;if\;\beta>0\\\
C(t+1)\;if\;\beta=0\\\ C\;if\;\beta<0.\end{array}\right.$

See Hoshino and Yamada [6], for instance, for the proof of Lemma 4.

######  Lemma 5.

Let $\mu,$ $\nu,$ $\tau$ and $z>0$, then

$z^{1-\nu}\int\limits_{0}^{z}(z-\xi)^{\nu-1}\xi^{\mu-1}e^{-\tau\xi}d\xi\leq
C\tau^{-\mu}$

where $C$ is a constant independent of $z.$

See Michalski [13] or Kirane and Tatar [8] for the proof of Lemma 5.

Let $y:[a,b]\to[0,\infty)$ be a continuous function satisfying

$y(t)\leq M+\int\limits_{a}^{t}\,\lambda(s)\,g(y(s))ds,\;a\leq t\leq b$ (2.1)

where $M>0$, $\lambda:[a,b]\to[0,\infty)$ is continuous and
$g:[0,\infty)\to[0,\infty)$ is continuous nondecreasing with

$g(y)>0,\quad\forall\;\;\;y>0.$

We define

$G_{M}(y)=\int\limits_{M}^{y}\frac{ds}{g(s)},\;y>0.$ (2.2)

The following Bihari-type inequality is classical and can be found in many
references. See for instance [16].

######  Lemma 6.

Under the above assumptions on $y,\lambda,g$, we have

$y(t)\leq G_{M}^{-1}\bigg{(}\int\limits_{a}^{t}\,\lambda(s)\,ds\bigg{)},$
(2.3)

for $a\leq t\leq T\leq b$ with $T>a$ satisfies

$\int\limits_{a}^{T}\,\lambda(s)\,ds<\int\limits_{M}^{\infty}\,\frac{ds}{g(s)}\,ds.$
(2.4)

######  Remark 2.

The classical Gronwall’s inequality corresponds to $g(y)=y$.

##  3\. Proof of Theorem 1

Define $\mathbf{z}=(u,v,w)$, we rewrite (1.6) as

$\displaystyle\left\\{\begin{array}[]{ccl}u_{t}-d_{1}\Delta\,u&=&f_{1}(t,\mathbf{z}),\\\
v_{t}-d_{2}\Delta v&=&f_{2}(t,\mathbf{z}),\\\
w_{t}-d_{3}\Delta\,w&=&f_{3}(t,\mathbf{z}),\\\ y_{t}-\Delta
y&=&f_{4}(t,\mathbf{z}),\;\;y(0)=0,\end{array}\right.$ (3.5)

where $f_{4}:=-f_{1}-f_{2}-f_{3}$ and the last equation in $y$ is added to
fulfill technical requirements as it will be clear below. From (3.5) we deduce
that

$\left(u+v+w+y\right)_{t}-\Delta\left(d_{1}u+d_{2}v+d_{3}w+y\right)=0.$ (3.6)

It follows that

$u(x,t)+v(x,t)+w(x,t)+y(x,t)-\Delta\mathbf{m}(x,t)=g(x),$ (3.7)

where

$\mathbf{m}(x,t)=\int_{0}^{t}\,\left[d_{1}u(x,s)+d_{2}v(x,s)+d_{3}w(x,s)+y(x,s)\right]\,ds,$
(3.8)

and

$g(x)=u_{0}(x)+v_{0}(x)+w_{0}(x).$ (3.9)

Therefore

$\|\mathbf{z}\|_{L^{\infty}_{x,t}}\leq\|\Delta\mathbf{m}\|_{L^{\infty}_{x,t}}+\|g\|_{L^{\infty}_{x}}\quad\text{and}\quad\|\mathbf{m}\|_{L^{\infty}_{x,t}}\leq\mathcal{C}_{T},$
(3.10)

where $\mathcal{C}_{T}$ depends continuously in $T>0$. Using (3.7), one can
write

$B(x,t)\partial_{t}\mathbf{m}-\Delta\mathbf{m}=g(x),$ (3.11)

where

$0<\underline{B}\leq
B(x,t)=\frac{u(x,t)+v(x,t)+w(x,t)+y(x,t)}{d_{1}u(x,t)+d_{2}v(x,t)+d_{3}w(x,t)+y(x,t)}\leq\overline{B}<\infty.$

Lemma 3.5 in [2] implies

$|\mathbf{m}(x,t)-\mathbf{m}(x^{\prime},t)|\leq\mathcal{C}_{T}|x-x^{\prime}|^{\delta},\;\;\;\forall\;0<t<T,$
(3.12)

for some $\delta\in(0,1)$. Equation (3.7) also implies

$\partial_{t}\mathbf{m}-\Delta\mathbf{m}=(d_{1}-1)u+(d_{2}-1)v+(d_{3}-1)w+g(x).$
(3.13)

By using the Hölder continuity of $\mathbf{m}$ and [2, Lemma 1.1, p. 283], we
infer

$|\nabla\mathbf{m}|\leq\mathcal{C}_{T}\,|\mathbf{z}|^{\frac{1-\delta}{2-\delta}}.$
(3.14)

Applying [2, Lemma 1.1, p. 283] to each equation of $q\in\\{u,v,w,y\\}$ leads
to

$|\nabla\,q|\leq\mathcal{C}_{T}\,\left(1+|\mathbf{z}|\right)^{\frac{3}{2}}.$
(3.15)

Using [2, Lemma 3.9, p. 295] together with (3.14) and (3.15) gives

$\displaystyle|\Delta\mathbf{m}|$ $\displaystyle\leq$
$\displaystyle\mathcal{C}_{T}\,|\nabla\mathbf{m}|^{1/2}\,\left(|\nabla\,u|+|\nabla\,v|+|\nabla\,w|+\,|\nabla\,g|\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\mathcal{C}_{T}\,|\mathbf{z}|^{\frac{1-\delta}{2(2-\delta)}}\,\left(1+|\mathbf{z}|^{3/2}\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\mathcal{C}_{T}\,\left(1+|\mathbf{z}|^{\frac{1-\delta}{2(2-\delta)}+\frac{3}{4}}\right).$

From the above estimate and (3.10), we obtain that

$\|\mathbf{z}\|_{L^{\infty}_{x,t}}\leq\mathcal{C}_{T}\,\left(1+\|\mathbf{z}\|_{L^{\infty}_{x,t}}^{\frac{1-\delta}{2(2-\delta)}+\frac{3}{4}}\right).$

Since $\frac{1-\delta}{2(2-\delta)}+\frac{3}{4}<1$ then
$\|\mathbf{z}\|_{L^{\infty}_{x,t}}\leq\mathcal{C}_{T}$. Consequently
$T_{max}=\infty$.

##  4\. Proof of Theorem 3

###### Proof.

A. The decay rate of $\left\|u\right\|_{C^{\mu}(\bar{\Omega})}$ and
$\left\|w\right\|_{C^{\mu}(\bar{\Omega})}$:

Clearly we have the integral equations associated with the first and third
equation of (1.6)

$u(t)=e^{-tA_{p}}u(0)+\int\limits_{0}^{t}e^{-(t-s)A_{p}}\left\\{a_{1}(s)w^{m}-a_{2}(s)u^{n}v^{k}\right\\}ds$
(4.1) $\begin{array}[]{c}w(t)=e^{-tG_{p}}w(0)\\\
+\int\limits_{0}^{t}e^{-(t-s)G_{p}}\left\\{-(a_{1}(s)+a_{3}(s))w^{m}-a_{4}(s)w+a_{2}(s)u^{n}v^{k}\right\\}ds.\end{array}$
(4.2)

Applying $A_{p}^{\alpha}$ and $G_{p}^{\alpha}$, $0<\alpha<1$ to both sides of
(4.1) and (4.2) respectively, we infer from Lemma 5 that

$\begin{array}[]{c}\left\|A_{p}^{\alpha}u\right\|_{p}\leq
C_{1}t^{-\alpha}e^{-bt}\left\|u_{0}\right\|_{p}+C_{1}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-b(t-s)}s^{\sigma_{1}}h_{1}(s)\left\|w^{m}\right\|_{p}ds\\\
+C_{1}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-b(t-s)}s^{\sigma_{2}}h_{2}(s)\left\|u^{n}v^{k}\right\|_{p}ds\end{array}$

and

$\begin{array}[]{ll}\left\|G_{p}^{\alpha}w\right\|_{p}\leq&C_{3}t^{-\alpha}e^{-bt}\left\|w_{0}\right\|_{p}\\\
&+C_{3}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-b(t-s)}(s^{\sigma_{1}}h_{1}(s)+s^{\sigma_{3}}h_{3}(s))\left\|w^{m}\right\|_{p}ds\\\
&+C_{3}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-b(t-s)}h_{4}(s)\left\|w\right\|_{p}ds\\\
&+C_{3}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-b(t-s)}a_{2}(s)\left\|u^{n}v^{k}\right\|_{p}ds.\end{array}$

As $0\leq v\leq M,$ for some $M>0,$ we have for all $t\geq\delta>0$

$\begin{array}[]{c}e^{bt}\left\|A_{p}^{\alpha}u\right\|_{p}\leq
C_{1}\delta^{-\alpha}\left\|u_{0}\right\|_{p}+C_{1}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{bs}s^{\sigma_{1}}h_{1}(s)\left\|w^{m}\right\|_{p}ds\\\
+C_{1}M^{k}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{bs}s^{\sigma_{2}}h_{2}(s)\left\|u^{n}\right\|_{p}ds,\end{array}$
(4.3)

and

$\begin{array}[]{c}e^{bt}\left\|G_{p}^{\alpha}w\right\|_{p}\leq
C_{3}\delta^{-\alpha}\left\|w_{0}\right\|_{p}+C_{3}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{bs}h_{4}(s)\left\|w\right\|_{p}ds\\\
+C_{3}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{bs}(s^{\sigma_{1}}h_{1}+s^{\sigma_{3}}h_{3})(s)\left\|w^{m}\right\|_{p}ds\\\
+C_{3}M^{k}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{bs}s^{\sigma_{2}}h_{2}(s)\left\|u^{n}\right\|_{p}ds.\end{array}$
(4.4)

Lemma 6 and the uniform boundedness of $w$ allow us to write

$\left\|w^{m}\right\|_{p}=\left\|w\right\|_{mp}^{m}\leq
C\left\|G_{p}^{\alpha}w\right\|_{p}^{m\theta}.\left\|w\right\|_{p}^{m(1-\theta)}\leq
C_{4}\left\|G_{p}^{\alpha}w\right\|_{p}^{m\theta}$ (4.5)

with $\frac{N(m-1)}{2pm\alpha}<\theta<1.$

In the same way we have

$\left\|u^{n}\right\|_{p}\leq
C_{5}\left\|A_{p}^{\alpha}u\right\|_{p}^{n\theta}\;with\;\frac{N(n-1)}{2pn\alpha}<\theta<1.$
(4.6)

Let us choose $l$ such that

$\max\left\\{1,\frac{N(m-1)}{2p\alpha},\frac{N(n-1)}{2p\alpha}\right\\}<l<\min\\{m,n\\}$

and $\theta=\frac{l}{m}$ in (4.5) and $\theta=\frac{l}{n}$ in (4.6). Then

$\left\|w^{m}\right\|_{p}\leq
C_{4}\left\|G_{p}^{\alpha}w\right\|_{p}^{l}\;and\;\left\|u^{n}\right\|_{p}\leq
C_{5}\left\|A_{p}^{\alpha}u\right\|_{p}^{l}.$ (4.7)

Taking (4.7) into account in (4.3) and (4.4), we obtain

$\begin{array}[]{c}e^{bt}\left\|A_{p}^{\alpha}u\right\|_{p}\leq
C_{1}\delta^{-\alpha}\left\|u_{0}\right\|_{p}\\\
+C_{6}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{b(1-l)s}s^{\sigma_{1}}h_{1}(s)\left(e^{bs}\left\|G_{p}^{\alpha}w\right\|_{p}\right)^{l}ds\\\
+C_{7}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{b(1-l)s}s^{\sigma_{2}}h_{2}(s)\left(e^{bs}\left\|A_{p}^{\alpha}u\right\|_{p}\right)^{l}ds.\end{array}$
(4.8)

and

$\begin{array}[]{c}e^{bt}\left\|G_{p}^{\alpha}w\right\|_{p}\leq
C_{3}\delta^{-\alpha}\left\|w_{0}\right\|_{p}+C_{3}\int\limits_{0}^{t}(t-s)^{-\alpha}h_{4}(s)e^{bs}\left\|w\right\|_{p}ds\\\
+C_{8}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{b(1-l)s}(s^{\sigma_{1}}h_{1}+s^{\sigma_{3}}h_{3})(s)\left(e^{bs}\left\|G_{p}^{\alpha}w\right\|_{p}\right)^{l}ds\\\
+C_{9}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{b(1-l)s}s^{\sigma_{2}}h_{2}(s)\left(e^{bs}\left\|A_{p}^{\alpha}u\right\|_{p}\right)^{l}ds.\end{array}$
(4.9)

The second term in the right hand side of (4.9) may be estimated in the
following manner, for $0<\varepsilon<b$

$\begin{array}[]{c}\int\limits_{0}^{t}(t-s)^{-\alpha}h_{4}(s)e^{bs}\left\|w\right\|_{p}ds=\int\limits_{0}^{t}(t-s)^{-\alpha}e^{\varepsilon
s}e^{(b-\varepsilon)s}h_{4}(s)\left\|w\right\|_{p}ds\\\
\leq\left(\int\limits_{0}^{t}(t-s)^{-\alpha q}e^{\varepsilon
qs}ds\right)^{\frac{1}{q}}\left(\int\limits_{0}^{t}h_{4}^{q^{*}}(s)\left(e^{(b-\varepsilon)s}\left\|w\right\|_{p}\right)^{q^{*}}ds\right)^{\frac{1}{q^{*}}}\\\
\leq C_{10}e^{\varepsilon
t}\left(\int\limits_{0}^{t}h_{4}^{q^{*}}(s)\left(e^{(b-\varepsilon)s}\left\|w\right\|_{p}\right)^{q^{*}}ds\right)^{\frac{1}{q^{*}}}.\end{array}$
(4.10)

We have used the Hölder inequality, Lemma 7 and the embedding
$D(G_{p}^{\alpha})\subset L^{p}$ to derive the last inequalities in (4.10).
Multiplying (4.8) and (4.9) by $e^{-\varepsilon t},$ setting

$U(t)=e^{(b-\varepsilon)t}\left\|A_{p}^{\alpha}u\right\|_{p},\text{
}W(t)=e^{(b-\varepsilon)t}\left\|G_{p}^{\alpha}w\right\|_{p}$

and taking into account (4.10) we find for $t\geq\delta>0$

$\begin{array}[]{c}U(t)\leq
C_{1}\delta^{-\alpha}\left\|u_{0}\right\|_{p}+C_{6}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-(b-\varepsilon)(l-1)s}s^{\sigma_{1}}h_{1}(s)W(s)^{l}ds\\\
+C_{7}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-(b-\varepsilon)(l-1)s}s^{\sigma_{2}}h_{2}(s)U(s)^{l}ds,\end{array}$

and

$\begin{array}[]{c}W(t)\leq
C_{3}\delta^{-\alpha}\left\|w_{0}\right\|_{p}+C_{10}\left(\int\limits_{0}^{t}h_{4}^{q^{*}}(s)W(s)^{q^{*}}ds\right)^{\frac{1}{q^{*}}}\\\
+C_{8}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-(b-\varepsilon)(l-1)s}(s^{\sigma_{1}}h_{1}+s^{\sigma_{3}}h_{3})(s)W(s)^{l}ds\\\
+C_{9}\int\limits_{0}^{t}(t-s)^{-\alpha}e^{-(b-\varepsilon)(l-1)s}s^{\sigma_{2}}h_{2}(s)U(s)^{l}ds.\end{array}$

Next, as $1-q\alpha>0,$ then using the Hölder inequality and Lemma 8, we
obtain for $t\geq\delta>0$

$\begin{array}[]{c}U(t)\leq
C_{1}\delta^{-\alpha}\left\|u_{0}\right\|_{p}+C_{11}\delta^{-\alpha}\left(\int\limits_{0}^{t}h_{1}^{q^{*}}(s)W(s)^{lq^{*}}ds\right)^{\frac{1}{q^{*}}}\\\
+C_{12}\delta^{-\alpha}\left(\int\limits_{0}^{t}h_{2}^{q^{*}}(s)U(s)^{lq^{*}}ds\right)^{\frac{1}{q^{*}}}\end{array}$

and

$\begin{array}[]{c}W(t)\leq
C_{3}\delta^{-\alpha}\left\|w_{0}\right\|_{p}+C_{10}\left(\int\limits_{0}^{t}h_{4}^{q^{*}}(s)W(s)^{q^{*}}ds\right)^{\frac{1}{q^{*}}}\\\
+C_{13}\delta^{-\alpha}\left(\int\limits_{0}^{t}h_{1}^{q^{*}}(s)W(s)^{lq^{*}}ds\right)^{\frac{1}{q^{*}}}+C_{14}\delta^{-\alpha}\left(\int\limits_{0}^{t}h_{3}^{q^{*}}(s)W(s)^{lq^{*}}ds\right)^{\frac{1}{q^{*}}}\\\
+C_{15}\delta^{-\alpha}\left(\int\limits_{0}^{t}h_{2}^{q^{*}}(s)U(s)^{lq^{*}}ds\right)^{\frac{1}{q^{*}}}.\end{array}$

Now we use the inequality

$(x_{1}+x_{2}+...+x_{n})^{r}\leq n^{r-1}(x_{1}^{r}+x_{2}^{r}+...+x_{n}^{r})$

to get

$\begin{array}[]{c}U(t)^{q^{*}}\leq
3^{q^{*}-1}\left(C_{1}\delta^{-\alpha}\left\|u_{0}\right\|_{p}\right)^{q^{*}}+3^{q^{*}-1}\left(C_{11}\delta^{-\alpha}\right)^{q^{*}}\int\limits_{0}^{t}h_{1}^{q^{*}}(s)W(s)^{lq^{*}}ds\\\
+3^{q^{*}-1}\left(C_{12}\delta^{-\alpha}\right)^{q^{*}}\int\limits_{0}^{t}h_{2}^{q^{*}}(s)U(s)^{lq^{*}}ds,\end{array}$
(4.11)

and

$\begin{array}[]{c}W(t)^{q^{*}}\leq
5^{q^{*}-1}\left(C_{3}\delta^{-\alpha}\left\|w_{0}\right\|_{p}\right)^{q^{*}}+5^{q^{*}-1}C_{10}^{q^{*}}\int\limits_{0}^{t}h_{4}^{q^{*}}(s)W(s)^{q^{*}}ds\\\
+5^{q^{*}-1}\left(C_{13}\delta^{-\alpha}\right)^{q^{*}}\int\limits_{0}^{t}h_{1}^{q^{*}}(s)W(s)^{lq^{*}}ds+5^{q^{*}-1}\left(C_{14}\delta^{-\alpha}\right)^{q^{*}}\int\limits_{0}^{t}h_{3}^{q^{*}}(s)W(s)^{lq^{*}}ds\\\
+5^{q^{*}-1}\left(C_{15}\delta^{-\alpha}\right)^{q^{*}}\int\limits_{0}^{t}h_{2}^{q^{*}}(s)U(s)^{lq^{*}}ds.\end{array}$
(4.12)

Putting $F(t)=U(t)^{q^{*}}+W(t)^{q^{*}},$ we infer from (4.11) and (4.12) that

$F(t)\leq
C_{0}\left(\left\|u_{0}\right\|_{p},\left\|w_{0}\right\|_{p}\right)+\int\limits_{0}^{t}h(s)(F(s)+F(s)^{l})ds,\;t\geq\delta>0$

where
$C_{0}\left(\left\|u_{0}\right\|_{p},\left\|w_{0}\right\|_{p}\right)=3^{q^{*}-1}\left(C_{1}\delta^{-\alpha}\left\|u_{0}\right\|_{p}\right)^{q^{*}}+5^{q^{*}-1}\left(C_{3}\delta^{-\alpha}\left\|w_{0}\right\|_{p}\right)^{q^{*}}$
and

$h(s)=\max\left\\{(C_{17}+C_{19})h_{1}^{q^{*}}(s),C_{20}h_{3}^{q^{*}}(s),(C_{18}+C_{21})h_{2}^{q^{*}}(s),C_{22}h_{4}^{q^{*}}(s)\right\\}$

where $C_{i}$, $i=17,18,...,22$ are the coefficients of the integral terms in
(4.11) and (4.12) in the order.

Let $G(z)=\int\limits_{z_{0}}^{z}\frac{dy}{y+y^{l}}.$ Then, by Lemma 6 we may
conclude that

$\begin{array}[]{c}F(t)\leq
G^{-1}\left[G\left(C_{0}\left(\left\|u_{0}\right\|_{p},\left\|w_{0}\right\|_{p}\right)\right)+\int\limits_{0}^{t}h(s)ds\right]\\\
\leq
C_{0}\left(1+C_{0}^{l-1}\right)^{\frac{1}{1-l}}e^{\int\limits_{0}^{t}h(s)ds}.\left[1-\frac{C_{0}}{1+C_{0}^{l-1}}e^{(l-1)\int\limits_{0}^{t}h(s)ds}\right]^{\frac{1}{1-l}}.\end{array}$

From our assumptions on $\left\|u_{0}\right\|_{p},$ $\left\|w_{0}\right\|_{p}$
and $h(t)$ we deduce that

$\left\|A_{p}^{\alpha}u\right\|_{p}\leq
Ce^{-(b-\varepsilon)t}\left(\left\|u_{0}\right\|_{p}+\left\|w_{0}\right\|_{p}\right)$
(4.13)

and

$\left\|G_{p}^{\alpha}w\right\|_{p}\leq
Ce^{-(b-\varepsilon)t}\left(\left\|u_{0}\right\|_{p}+\left\|w_{0}\right\|_{p}\right).$
(4.14)

The decay rates in $C^{\mu}(\bar{\Omega}),$ $\mu\in[0,2)$ follow from Lemma 4.

B. The decay rate of $\left\|v-v_{\infty}\right\|_{C^{\mu}(\bar{\Omega})}:$

As in Hoshino [5], let us write

$v-v_{\infty}=\left(Q_{0}v(t)-v_{\infty}\right)+Q_{+}v(t),$

and estimate the terms in the right hand side separately.

a. The estimation of $Q_{0}v(t)-v_{\infty}:$

Integrating the second equation in (1.6) over $(0,t)\times\Omega$, we have

$\int\limits_{\Omega}v(x,t)dx+\int\limits_{0}^{t}\int\limits_{\Omega}\left\\{(a_{1}(s)+a_{3}(s))w^{m}-a_{2}(s)u^{n}v^{k}\right\\}dxds=\int\limits_{\Omega}v_{0}(x)dx.$

It appears then that

$\left|Q_{0}v(t)-v_{\infty}\right|=\left|\Omega\right|^{-1}\left|\int\limits_{t}^{\infty}\int\limits_{\Omega}\left\\{(a_{1}(s)+a_{3}(s))w^{m}-a_{2}(s)u^{n}v^{k}\right\\}dxds\right|.$

In the rest of the proof, $C$ will denote a generic positive constant which
may be different at different occurrences.

Using (4.7), (4.13), (4.14) and the Hölder inequality we see that

$\begin{array}[]{l}\left|Q_{0}v(t)-v_{\infty}\right|\\\ \leq
C\int\limits_{t}^{\infty}\left\\{(s^{\sigma_{1}}h_{1}+s^{\sigma_{3}}h_{3})(s)e^{-l(b-\varepsilon)s}+s^{\sigma_{2}}h_{2}(s)M^{k}e^{-l(b-\varepsilon)s}\right\\}ds\\\
\leq
C\int\limits_{t}^{\infty}\left\\{\sum\limits_{i=1}^{3}s^{\sigma_{i}}h_{i}(s)\right\\}e^{-l(b-\varepsilon)s}ds\\\
\leq
C\int\limits_{t}^{\infty}\left\\{\sum\limits_{i=1}^{3}s^{\sigma_{i}}h_{i}(s)\right\\}e^{-(b-\varepsilon)s}.e^{-(l-1)(b-\varepsilon)s}ds\\\
\leq
Ce^{-(b-\varepsilon)t}\int\limits_{t}^{\infty}\left\\{\sum\limits_{i=1}^{3}s^{\sigma_{i}}h_{i}(s)\right\\}e^{-(l-1)(b-\varepsilon)s}ds\\\
\leq Ce^{-(b-\varepsilon)t}.\end{array}$ (4.15)

We also used Lemma 7 in the last inequality.

b. The estimation of $Q_{+}v(t):$

In order to estimate $Q_{+}v(t)$ let us apply $B_{p+}^{\alpha}Q_{+}$ to the
integral equation associated with the second equation in (1.6). We find for
all $t\geq\delta>0$

$\begin{array}[]{c}B_{p+}^{\alpha}Q_{+}v(t)=B_{p+}^{\alpha}e^{-(t-\delta)B_{p+}}Q_{+}v(\delta)\\\
+\int\limits_{\delta}^{t}B_{p+}^{\alpha}e^{-(t-s)B_{p+}}Q_{+}\left[(a_{1}(s)+a_{3}(s))w^{m}-a_{2}(s)u^{n}v^{k}\right]ds.\end{array}$

Taking the $L^{p}$-norm and using the second inequality in Lemma 5, we obtain

$\begin{array}[]{c}\left\|B_{p+}^{\alpha}Q_{+}v(t)\right\|_{p}\leq
C_{2}(t-s)^{-\alpha}e^{-d_{2}(t-s)\lambda}\left\|Q_{+}v(\delta)\right\|_{p}+C_{2}\left\|Q_{+}\right\|_{L^{p}(\Omega)\rightarrow
L^{p}(\Omega)}\\\
\int\limits_{\delta}^{t}(t-s)^{-\alpha}e^{-d_{2}\lambda(t-s)}\left\\{(s^{\sigma_{1}}h_{1}+s^{\sigma_{3}}h_{3})\left\|w^{m}\right\|+s^{\sigma_{2}}h_{2}M^{k}\left\|u^{n}\right\|\right\\}ds.\end{array}$

Next, using (4.7), (4.13) and (4.14) we see that for all $t\geq\delta+T$

$\begin{array}[]{ll}\left\|B_{p+}^{\alpha}Q_{+}v(t)\right\|_{p}&\leq
Ce^{-d_{2}\lambda t}\left\|Q_{+}v(\delta)\right\|_{p}\\\ &+Ce^{-d_{2}\lambda
t}\int\limits_{\delta}^{t}(t-s)^{-\alpha}e^{d_{2}\lambda
s}e^{-l(b-\varepsilon)s}\left(\sum\limits_{i=1}^{3}s^{\sigma_{i}}h_{i}(s)\right)ds\\\
&\leq Ce^{-d_{2}\lambda
t}\left\\{\left\|Q_{+}v(\delta)\right\|_{p}+\int\limits_{0}^{t-\delta}(t-\delta-s)^{-\alpha}\right.\\\
&\left.\times
e^{-\left[l(b-\varepsilon)-d_{2}\lambda\right](s+\delta)}\left(\sum\limits_{i=1}^{3}(s+\delta)^{\sigma_{i}}h_{i}(s+\delta)\right)ds\right\\}.\end{array}$
(4.16)

  * (i)

If $d_{2}\lambda<lb$, then choose $\varepsilon$ such that $0<l\varepsilon<lb-
d_{2}\lambda.$ Hence, we may apply Lemma 8, together with the Hölder
inequality to get

$\left\|B_{p+}^{\alpha}Q_{+}v(t)\right\|_{p}\leq Ce^{-d_{2}\lambda
t}\;for\;all\;t\geq\delta+T.$ (4.17)

  * (ii)

If $d_{2}\lambda\geq lb,$ then $l(b-\varepsilon)-d_{2}\lambda<0.$ Multiplying
by $e^{\rho(s+\delta)}.e^{-\rho(s+\delta)},$ with
$\rho>d_{2}\lambda-l(b-\varepsilon),$ the integrand in (4.16), using the
Hölder inequality and Lemma 8, we see that

$\displaystyle\left\|B_{p+}^{\alpha}Q_{+}v(t)\right\|_{p}$ $\displaystyle\leq$
$\displaystyle Ce^{-d_{2}\lambda
t}\left\\{\left\|Q_{+}v(\delta)\right\|_{p}+e^{-\tilde{\rho}t}\right\\}$
$\displaystyle\leq$ $\displaystyle Ce^{-(d_{2}\lambda-\tilde{\rho})t}.$

The conclusion follows from (4.15), (4.17) and ((ii)).

∎

##  5\. Concluding Remarks

Some remarks are in order:

  * •

The condition $a_{1}(t)\leq Ca_{3}(t)$ needed in Theorem 2 has not been used
in Theorem 3. So, the decay rates hold provided one may prove a convergence
result without this condition.

  * •

The assumption
$\max\left\\{1,\frac{N(m-1)}{2p\alpha},\frac{N(n-1)}{2p\alpha}\right\\}<\min\\{m,n\\}$
may be relaxed somewhat using different $l_{1}$ and $l_{2}$ and applying Lemma
9 with $p=2.$

  * •

We also have exponential decay in (b) (ii) without the growth condition on
$h_{i}$, $i=1,2,3$ in case $\sigma_{i}=0,$ $i=1,2,3.$

  * •

It is possible to obtain sharper estimates using comparison results by
replacing the bounds $M^{k}$ with $v_{\infty}^{k}-Ce^{-(b-\varepsilon)t}$ in
(4.3) and $v_{\infty}^{k}+Ce^{-(b-\varepsilon)t}$ in (4.4), for large values
of $t$.

## Acknowledgment

We express our deep thanks to Bao Quoc Tang for giving us the proof of Theorem
1 and for interesting discussions. The second author would like to thank King
Fahd University of Petroleum and Minerals (Saudi Arabia) for its support.

## References

  * [1] H. Amann, Global existence for semilinear parabolic systems, J. Reine Angew. Math., 360 (1985), 47–83.
  * [2] K. Fellner, J. Morgan and Bao Q. Tang, Global classical solutions to quadratic systems with mass control in arbitrary dimensions, Ann. Inst. Henri Poincaré, Anal. Non Linéaire, 37 (2020), 281–307.
  * [3] W. Feng, Coupled system of reaction-diffusion equations and applications in Carrier faciliated diffusion, Nonlinear Anal. T.M.A., Vol. 17, No. 3 (1991), 285-311.
  * [4] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, Vol. 840, Springer-Verlag, Berlin-New York, 1981.
  * [5] H. Hoshino, On the convergence properties of global solutions for some reaction diffusion systems under Neumann boundary conditions, Diff. and Integ. Eqs., Vol. 9, No. 4 (1996), 761-778.
  * [6] H. Hoshino and Y. Yamada, Solvability and smoothing effect for semilinear parabolic equations, Funkc. Ekvac., 34 (1994), 475-494.
  * [7] C. S. Kahane, On the asymptotic behavior of solutions of nonlinear parabolic systems under Robin type boundary conditions, Funkc. Ekvac., 26 (1983), 51-78.
  * [8] M. Kirane and N.-E. Tatar, Global existence and stability of some semilinear problems, Arch. Math. Brno, Tomus 36 (2000), 1-12.
  * [9] M. Kirane and N.-E. Tatar, Convergence rates for a reaction diffusion system, Zeit. Anal. Anw. (J. Anal.Math.), 20 (2001) 2, 347-357.
  * [10] E. H. Laamri and M. Pierre, Global existence for reaction-diffusion systems with nonlinear diffusion and control of mass, Ann. Inst. Henri Poincaré, Anal. Non Linéaire, 34 (2017), 571–591.
  * [11] E. H. Laamri and B. Perthame, Reaction-diffusion systems with initial data of low regularity, J. Differential Equations, 269 (2020), 9310–9335.
  * [12] R. H. Martin and M. Pierre, Nonlinear reaction-diffusion systems, in Nonlinear Equations in the Applied Sciences, Math. Sci. Engrg., 185 W. F. Ames and C. Rogers, Editors, Academic Press, Boston, MA 1992, pp. 363-398.
  * [13] M. W. Michalski, Derivatives of non integer order and their applications, Dissertationes Mathematicae, Polska Akademia Nauk, Instytut Matematyczny, Warszawa, (1993).
  * [14] J. Morgan, Global existence for semilinear parabolic systems, SIAM J. Math. Anal., Vol. 20, No. 5 (1989), 1128-1144.
  * [15] M. Pierre and D. Schmitt, Blowup in reaction-diffusion systems with dissipation of mass, SIAM J. Math. Anal., 28 (1997), 259–269.
  * [16] M. Pinto, Integral inequalities of Bihari-type and applications, Funkc. Ekvac., 33, No. 3 (1990), 387-404.
  * [17] F. Rothe, Global Solutions of Reaction-Diffusion Systems, Lecture Notes in Mathematics, Vol. 1072, Springer-Verlag, Berlin, 1984.
  * [18] W. Ruan, Bounded solutions for reaction diffusion systems with nonlinear boundary conditions, Nonlinear Anal. T.M.A., Vol. 14, No. 12 (1990), 1051-1077.
  * [19] M. X. Wang, Asymptotic behavior of solutions to some reaction diffusion systems, Chinese J. Contemp. Math., Vol. 18, No. 3 (1997), 249-260.


# Linear correlation between active and resistive stresses informs on force
generation and stress transmission in adherent cells

Hélène Delanoë-Ayari<EMAIL_ADDRESS>Univ. Claude Bernard
Lyon1, CNRS, Institut Lumière Matière, 69622 Villeurbanne, France  Nicolas
Bouchonville  Univ. Grenoble Alps, CNRS, LTM, 38000 Grenoble, France  Marie
Courçon  Univ. Grenoble Alps, CEA, Inserm, BIG-BGE, 38000 Grenoble, France
Alice Nicolas<EMAIL_ADDRESS>Univ. Grenoble Alps, CNRS, LTM, 38000
Grenoble, France

###### Abstract

Animal cells are active, contractile objects. While bioassays address the
molecular characterization of cell contractility, the mechanical
characterization of the active forces in cells remains challenging. Here by
confronting theoretical analysis and experiments, we calculated both the
resistive and the active components of the intracellular stresses that build
up following cell adhesion. We obtained a linear relationship between the
divergence of the resistive stress and the traction forces, which we show is
the consequence of the cell adhering and applying forces on the surface only
through very localized adhesion points (whose size is inferior to our best
resolution, of 400 nm). This entails that there is no measurable forces
outside of these active point sources, and also that the resistive and active
stresses inside cells are proportional.

biological physics, mechanobiology, traction force microscopy, intracellular
stress microscopy, elasticity

###### pacs:

10xxx

Animal cells have contractile capabilities that make cells tensed objects.
This contractility allows adherent cells to probe the mechanical properties of
their environment and adapt to them [1, 2, 3]. Dysfunction of cell
contractility is a hallmark of many pathologies, such as cancers, cardiac or
brain pathologies [4, 5]. As it is strictly regulated and adapts to external
physical or chemical perturbations [6], the analysis of cell contractility
often brings information on the interplay of specific signaling pathways with
the extracellular environment. For example, stem cell differentiation was
shown to be closely regulated by the level of contractility of the tissue they
are part of [7]. When asking about cell contractility, the biological question
is in general to identify, locate and quantify the biochemical processes in
cells that give rise to cellular forces, contractile or tensile stresses. The
activity of molecular motors for instance results in mechanical stresses [2,
8]. Changes in the conformation of these proteins generate molecular movements
that mechanically translate into generation of forces at the molecular level.
In cell biology, these sources of stress are sought using molecular markers
thus setting assumptions on the biological nature of the intracellular stress
generators.

Figure 1: Modeling of an adherent cell for intracellular stress calculation.
a) Schematic of an adherent cell. Acto-myosin filaments (in red) are attached
to focal adhesions (in green) and may raise tension in cell body. $h$ is the
thickness of the layer where the stresses transmitted to the substrate are
generated. b) Elastic model for a cell or a cell colony (in grey) firmly
adhered to a semi infinite deformable matrix (in blue). The respective Young’s
moduli for the thin film and the semi infinite layer are $E_{c}$ and $E_{m}$
and their Poisson’s ratio $\nu_{c}$ and $\nu_{m}$. The cell is assumed to bear
a point of stress generation, $\vec{f}_{act}$ (red square). The thin film
opposes a resistance $\vec{f}_{c}$ to the active stress, and the matrix
opposes $-\vec{f}_{m}$. c) Our results indicate that intracellular stresses
are transmitted to the substrate through discrete anchorages of size smaller
than the experimental sampling size $a$. This transmission could either come
from discrete connections of the stress generators (e.g. the acto-myosin
stress fibers and paxillin-stained adhesions) or of unstained adhesion sites
(depicted in yellow).

More recently, a need for label-free approaches to assess cell contractility
has emerged. Their objective is to identify the areas of stress generation and
to quantify their amplitude. Optical methods have been proposed that measure
the density of cytoskeleton fibers in the absence of staining [9]. With even
less assumption on the origin of stress generation, mechanical approaches have
been implemented that quantify intracellular mechanical stresses [10, 11, 12,
13]. These methods are based on the measurement of the deformation of the
extracellular environment the cells are adhering to and exploit it to
calculate cell internal stresses. Here we focus on these mechanical
approaches.

By combining them as described in [14], we observe a linear correlation
between the active and resistive components of the intracellular stress
tensors. Complementing this observation with theoretical approaches, we bring
a new picture of the interaction of the cells with the substrate, showing the
existence of discrete mechanical links between the stress generators and the
substrate at submicron scale.

There exist different techniques for calculating cellular stress in cells. One
set of methods is based on the writing of force conservation inside a 2D
material. It uses as input the traction forces $\vec{f}_{m}$ exerted by the
cells on its environment and solves the 2D equation:

$hdiv(S_{tot})=\vec{f}_{m}$ (1)

where $S_{tot}$ is the total intracellular stress, meaning, the sum of the
active and resistive stresses $S_{act}$ and $S_{c}$ that originate from the
active cellular surface forces $\vec{f}_{act}$ and the reaction surface forces
from the cell body, $\vec{f}_{c}$ (Fig. 1b). $h$ in Eq. (1) is the thickness
of the layer in the cell where the stresses transmitted to the substrate are
generated. Eq. (1) is solved by using either standard finite elements in the
way proposed by Tambe _et al._ [15] for Monolayer Stress Microscopy (MSM) or a
Bayesian approach as proposed by Nier _et al._ [13] for Bayesian Inference
Stress Microscopy (BISM). Both these techniques enable to recover the total
intracellular stress, $S_{tot}=S_{act}+S_{c}$. Differently, as shown in [14],
Intracellular Stress Microscopy (ISM) enables to recover the Young’s modulus-
normalized resistive stress tensor $S_{c}/E_{c}$. This latter technique does
not require the calculation of the surface forces $\vec{f}_{m}$ from the cells
to the substrate, but is based solely on the continuity of the displacement at
the cell/substrate interface [12]. For cell biology issues, a quantity of
prime interest is $\vec{f}_{act}$, the internal cellular surface forces at the
origin of cell contractility. In principle, combination of MSM or BISM and ISM
will provide $hS_{act}$ from which $\vec{f}_{act}$ can be derived [14]. We
thus decided to calculate both quantities $S_{tot}$ and $S_{c}$ using BISM and
ISM, and we present here an in-depth exploration of their relationship
obtained in two different cell types.

We first investigated the intracellular stresses in rat embryonic fibroblast
cell line REF52 (Fig. 2). The REF52 cell line we used was stably transfected
with fluorescent paxillin (gift from A. Bershadsky), so to compare the
location of intracellular stresses and paxillin-stained focal adhesions. The
geometry of the single cells was consistent with the plaque approximation, the
height of the cells being at a maximum of 5 $\mu$m (data not shown) to be
compared to their in-plane extent of order of 50 to 100 $\mu$m. Single cells
were grown on a soft polyacrylamide hydrogel of 3 kPa functionalized with
fibronectin. The hydrogel was loaded with a high density of 200 nm fluorescent
markers. The deformation field of the substrate was quantified by comparing
images of beads located close to its top surface in the presence of cells and
when the cells are removed. Beads displacements were measured using a
pyramidal optical flow algorithm (SI-1). The surface forces $\vec{f}_{m}$ were
calculated using fast Fourier transformation of the displacement field [16].
We first observed that the traction force field did not evidence correlations
with the distribution of the paxillin-stained adhesions (Fig. 2b). This
suggests that cell intracellular stresses are transmitted to the extracellular
matrix also out of these adhesions. This result is not surprising as Zamir et
al. have shown that in REF 52 cells paxillin staining does not stain tensin
rich adhesions [17]. This again promotes a label-free approach, as one can
never be sure that labeling one (or even several proteins) will guarantee the
observation of all sites of interest for active stress generation.

Since the calculation of $S_{c}$ by ISM only makes sense when the cell body is
firmly bound to the substrate, we limited stress calculation to paxillin-
positive regions and to regions where $\vec{f}_{m}$ is above the noise level
(Fig. 2c). In these regions, the fact that the traction force field
$\vec{f}_{m}$ is out of the noise implies that the cell is adhered and
intracellular stresses are transmitted to the substrate. $S_{tot}$ was
calculated in the same regions using BISM algorithm, following the methodology
described in [14] (Fig. 2d). Comparison of BISM and ISM revealed a linear
correlation between both, with a negative slope (Fig. 2e). In addition,
following a previous work where we had reported on a linear relationship
between the amplitudes of $divS_{c}$ and $\vec{f}_{m}$ [12], we confirmed this
linear correlation for this other cell type. Components of the divergence of
the resistive stress tensor $S_{c}$ correlate with surface force components
$\vec{f}_{m}$ with a minus sign (Fig. 2f):

$div\frac{S_{c}}{E_{c}}=-\frac{\vec{f}_{m}}{\ell E_{m}}$ (2)

with $E_{m}$ the Young’s modulus of the matrix and $\ell$ a characteristic
length. As visible on Fig. 2f, positively paxillin-labeled pixels are
indistinguishable from unlabeled pixels. This observation provides an
additional argument for enlarging the regions of cell adhesion out of
paxillin-positive adhesions [12].

Figure 2:  a) Focal adhesions in REF52 stably transfected for YFP-paxillin.
Bar 20 $\mu m$. b) Amplitude of $\vec{f}_{m}/E_{m}$ superimposed with the cell
contour and the contour of the paxillin-stained adhesions (in white) shows
significant stresses out of paxillin-stained adhesions. c) Amplitude of
$S_{c}/E_{c}$ measured at places where $\vec{f}_{m}$ exceeds noise level. d)
Amplitude of $hS_{tot}/E_{m}$ in $\mu m$, calculated at the same places
(regularization parameter $L=0.06$). e) The components of $hS_{tot}/E_{m}$ and
$S_{c}/E_{c}$ show a linear correlation (slope 2.28 $\mu m^{-1}$). f) The
components of $divS_{c}/E_{c}$ and $\vec{f}_{m}/E_{m}$ are proportional (slope
2.27 $\mu m^{-1}$). Green dots are for paxillin-labeled pixels, blue dots for
unlabeled pixels.

To understand these linear correlations, we calculated the theoretical
relationship between $divS_{c}$ and $\vec{f}_{m}$ in a model system that
consists of a thin elastic layer continuously bound to a semi-infinite elastic
medium and stressed by a local stress field (see Fig. 1b). As the surface
forces $\vec{f}_{m}$ are linked to the displacement field through a Green
function [18], the relationship between $divS_{c}$ and $\vec{f}_{m}$ is of
similar shape: a non local relationship, with a combined influence of the
stresses from both in-plane directions. We however obtained that this non
local relationship can be approximated to a local proportionality because (i)
the off-diagonal terms in the Green function are two orders of magnitude
smaller than the diagonal terms, and (ii) the diagonal terms are fast decaying
functions close to the force point (SI-2). Because of this fast decay, the
relationship between $divS_{c}$ and $\vec{f}_{m}$ is sensitive to the ratio of
the lateral extent of $\vec{f}_{m}$ and the sampling size of the grid that is
used to perform Traction Force Micoscopy (TFM) or stress calculations.
Actually, a linear correlation between $divS_{c}$ and $\vec{f}_{m}$ was only
obtained when the lateral extent of the surface forces $\vec{f}_{m}$ is
smaller than the sampling size (Fig. 3). The opposite case, where the
amplitude of the surface forces spreads on a width larger than the sampling
size leads to a non linear correlation, different from the experimental
observation (Fig. 3). Facing the model with the experimental observation thus
leads to the conclusion that the traction forces $\vec{f}_{m}$ apply on areas
that are smaller than the size of the sampling grid that is used in TFM. So
due to the size of the sampling, $\vec{f}_{m}$ appears as point forces. The
model then predicts:

$hdivS_{c}=\vec{f}_{c}\simeq-\alpha\vec{f}_{m}$ (3)

where $\alpha=\frac{\pi
hE_{c}(1+\nu_{m})(3-2\nu_{m}-\nu_{c})}{3aE_{m}(1-\nu_{c}^{2})}$ with $a$ the
size of the sampling grid and $\nu_{c}$ and $\nu_{m}$ the Poisson’s ratios of
the cell and the substrate (SI-2). For wider distributions of $\vec{f}_{m}$,
the correlation showed two branches (Fig. 3b), also observed in 3D FEM
simulation [14], a consequence of the oscillations of the Green function that
couples both quantities (Fig. S2). From this analysis, we could conclude that
the proportionality between $divS_{c}$ and $\vec{f}_{m}$ that we observe in
the experiment is indeed related to the small extent of the traction forces
compared to the sampling size, and is anyhow an approximate linearity.
Combined with the observation that the amplitude of $\vec{f}_{m}$ is above the
noise level in a large part of the cell (Fig. S5), we conclude that the
surface forces $\vec{f}_{m}$ are concentrated to very local areas whose size
is below our in-plane resolution of $0.7\mu m$, but are distributed almost
everywhere beneath the cell, not restricted to paxillin-stained adhesions.

Figure 3: Analysis of the correlation of $divS_{c}$ and $\vec{f}_{m}$. a)
Scheme of the protocol used for calculating curves drawn in (b). A dipole
force with a Gaussian distribution whose width $\sigma$ is either smaller
(resp. larger) than the sampling size $a$ (top, resp. bottom) is simulated for
calculating the divergence of Sc according to Eq. (S8). b) $divS_{c}$ and
$\vec{f}_{m}$ show a linear correlation when $\sigma<a$ ($a=8$ pixels, red:
$\sigma=1$ pixel; blue: $\sigma=25$ pixels). Goodness of the fit for the red
curve: $r^{2}=0.83$. The dark line is a bin average of the blue points.

Eq. (2) has introduced a characteristic length scale $\ell$ that should
compare to $h/\alpha$ in Eq. (3). We artificially reduced the in-plane
resolution to probe the dependency of $\ell$ with the sampling size $a$. As
shown in Fig. S3, we obtained that $\ell$ is proportional to $a$, as predicted
in Eq. (3) (SI-3). This confirmed our analysis on the role of the sampling
size in the relation between $divS_{c}$ and $\vec{f}_{m}$. This then implies
that at points where $\vec{f}_{m}\neq\vec{0}$,

$\vec{f}_{act}\simeq(1+\alpha)\vec{f}_{m}$ (4)

and

$S_{c}\simeq-\alpha S_{tot}+\Phi$ (5)

with $\Phi$ a zero divergence stress tensor set by the boundary conditions.
Both linear correlations Eqs (4) and (5) are a direct consequence of the
linear correlation between $divS_{c}$ and $\vec{f}_{m}$ (Eq. (3)) 111It should
be noted that linearity is optimal when the regularization parameter in BISM
calculation is chosen with the chi2 criterion (see SI-1 and Fig. S6). that we
also observed experimentally (Eq. (2)) 222In the experiment, the zero
divergence stress tensor $\Phi$ in Eq. (5) is not apparent. This could be that
the boundary conditions impose an amplitude for $\Phi$ that is buried in the
noise or that its contribution is filtered by the regularization step in the
calculation of $hS_{tot}$, that filters low frequencies (see SI-1).. Then,
altogether these results suggest that either the stress generators are small
entities whose size is smaller than our sampling size, or they are
mechanically linked to the substrate by discrete anchors whose size is smaller
than our sampling size, not restricted to focal adhesions (Fig. 1d). This
result is actually consistent with other studies that identified myosin or
acto-myosin microfilaments as stress generators with size of few hundreds of
nanometers, far below the present in-plane resolution [21, 22].

We wondered whether the linear correlation between $\vec{f}_{m}$ and
$divS_{c}$ we observed with REF 52 cells was specific to the experimental
conditions used here. To test the robustness of our observations, Dr.
Fritzsche’s group made available raw data obtained with Hela cells expressing
GFP-paxillin cultured on a 40 kPa polyacrylamide hydrogel loaded with
fluorescent beads of 40 nm diameter [23] (Fig. 4). The position of the beads
was imaged with STED microscopy, as described in [24]. In this experiment, the
pixel size is about 20 nm to be compared to 100 nm in our experiment. The much
stiffer substrate allowed cells to develop more mature focal adhesions,
although it may limit our capability to detect small stresses as small
deformations may be hidden by the noise. The enhanced resolution of STED
microscopy allowed us to reach a spatial resolution of 400 nm, 2000 beads
being successfully tracked in the image. Figure 4b shows the calculated
displacement field of the beads. Here, large traction forces were observed in
focal adhesions at the periphery of the cell (Fig. 4c). For the first time,
thanks to the enhanced resolution, alternating compressive and tensile
stresses were made visible within the focal adhesions (Fig. 4d), as predicted
by theories that address the growth of focal adhesions in the force direction
[25, 26]. As for the REF52 cells grown on a much softer substrate, a linear
correlation between $divS_{c}/Ec$ and $\vec{f}_{m}/E_{m}$ was observed (Fig.
4e). This confirms that this relation does not come from bias in the
experimental set-up. So we conclude that the measurement of the traction
forces $\vec{f}_{m}$ gives information on the location of the intracellular
stress generators (Fig. 1d).

Figure 4: a) Paxillin expressing Hela cells on 40 kPa polyacrylamide hydrogels
imaged with STED microscopy. Bar $5\mu m$. b) Displacement field obtained with
the KLT optical flow algorithm. The white arrow is $0.5\mu m$ long. c)
Traction forces concentrate in focal adhesions and show dotted patterns. White
lines delineate cell periphery and focal adhesions. d) Maps of the stress
components, superimposed with the contours of the cell and the focal adhesions
(white lines). Local compressive and tensile stresses are visible within focal
adhesions. e) The linear correlation between $\vec{f}_{m}/E_{m}$ and
$div(S_{c}/E_{c})$ is still observed with these mature focal adhesions. Raw
data are provided by H. Colin-York and M. Fritzsche.

In conclusion, we report on a linear correlation between the divergence of the
stress tensor in the cell body and the forces that are transmitted to the
substrate (Fig. 2f). We show that this linear correlation implies that the
transmission of the cellular stresses to the substrate is performed through
local links whose size is smaller than the sampling size of the experiment, as
depicted in Fig. 1d, and that this conclusion is independent on assumptions on
cell rheology (SI-4). When the cells have a linear elastic behavior, we show
that stress generation following cell adhesion leads as a first approximation
to the production of a proportional resistive stress in the cell body. Thus
quantification of the intracellular stresses either by MSM, BISM or ISM brings
similar qualitative results. It also makes it possible to localize stress
generators by measuring the surface forces $\vec{f}_{m}$ that cells transmit
to the extracellular environment thus highlighting the sensitivity and the
relevance of mechanical analysis as companion technique of biological
analysis.

###### Acknowledgements.

The authors are indebted to D. Gulino-Debrac for allowing them to use the
biology lab and to P. Marcq for the provision of the BISM calculation code. A.
N. and H. D.-A. deeply acknowledge H. Colin-York and M. Fritzsche for
providing additional raw data to test the accuracy of the analysis. This work
was initiated by very fruitful discussions with E. Mazza, L. Filotto, P.
Silberzan and T. Vourc’h. H. D. and A. N. are grateful to them. The authors
also thank F. Graner for critical reading. A. N. and N.B. acknowledge the
support by ANR-12-JSVE05-0008.

## References

  * Polte  _et al._ [2004] T. R. Polte, G. S. Eichler, N. Wang, and D. E. Ingber, Am. J. Physiol. Cell Physiol. 286, C518 (2004).
  * Doss  _et al._ [2020] B. L. Doss, M. Pan, M. Gupta, G. Grenci, R.-M. Mège, C. T. Lim, M. P. Sheetz, R. Voituriez, and B. Ladoux, Proc. Natl Acad. Sci. U.S.A. 117, 12817 (2020).
  * Zhang  _et al._ [2020] J. Zhang, F. Alisafaei, M. Nikolić, X. A. Nou, H. Kim, V. B. Shenoy, and G. Scarcelli, Small 16, 1907688 (2020).
  * Northcott  _et al._ [2018] J. M. Northcott, I. S. Dean, J. K. Mouw, and V. M. Weaver, Front. Cell Dev. Biol. 6, 17 (2018).
  * Javier-Torrent and Saura [2020] M. Javier-Torrent and C. A. Saura, Cells 9, 1926 (2020).
  * Bordeleau and Reinhart-King [2016] F. Bordeleau and C. A. Reinhart-King, F1000Research 5, 1819 (2016).
  * Ning  _et al._ [2021] W. Ning, A. Muroyama, H. Li, and T. Lechler, Cell Stem Cell 28, 436 (2021).
  * Peyton  _et al._ [2007] S. R. Peyton, C. M. Ghajar, C. B. Khatiwala, and A. J. Putnam, Cell Biochem. Biophys. 47, 300 (2007).
  * Wang  _et al._ [2018] W. Wang, J. P. Miller, S. C. Pannullo, C. A. Reinhart-King, and F. Bordeleau, J. Biophotonics 11, e201800008 (2018).
  * Wang  _et al._ [2002] N. Wang, I. M. Tolić-Nørrelykke, J. Chen, S. M. Mijailovich, J. P. Butler, J. J. Fredberg, and D. Stamenović, Am. J. Physiol. Cell Physiol. 282, 606 (2002).
  * Tambe  _et al._ [2011] D. T. Tambe, C. Corey Hardin, T. E. Angelini, K. Rajendran, C. Y. Park, X. Serra-Picamal, E. H. Zhou, M. H. Zaman, J. P. Butler, D. A. Weitz, J. J. Fredberg, and X. Trepat, Nat. Mat. 10, 469 (2011).
  * Moussus  _et al._ [2014] M. Moussus, C. der Loughian, D. Fuard, M. Courçon, D. Gulino-Debrac, H. Delanoë-Ayari, and A. Nicolas, Soft Matter 10, 2414 (2014).
  * Nier  _et al._ [2016] V. Nier, S. Jain, C. T. Lim, S. Ishihara, B. Ladoux, and P. Marcq, Biophys. J. 110, 1625 (2016).
  * Delanoë-Ayari and Nicolas [2022] H. Delanoë-Ayari and A. Nicolas, submitted to PRE (2022).
  * Tambe  _et al._ [2013] D. T. Tambe, U. Croutelle, X. Trepat, C. Y. Park, J. H. Kim, E. Millet, J. P. Butler, and J. J. Fredberg, PLoS One 8, e55172 (2013).
  * Butler  _et al._ [2002] J. P. Butler, I. Tolić-Nørrelykke, B. Fabry, and J. J. Fredberg, Am. J. Physiol. Cell Physiol. 282, C595 (2002).
  * Zamir  _et al._ [1999] E. Zamir, B. Z. Katz, S.-i. Aota, K. M. Yamada, B. Geiger, and Z. Kam, J. Cell Sci. 112, 1655 (1999).
  * Landau and Lifshitz [1970] L. D. Landau and E. M. Lifshitz,  _Theory of elasticity_ , 2nd ed., Course of Theoretical Physics, Vol. 7 (Pergamon Press, Oxford, 1970).
  * Note [1] It should be noted that linearity is optimal when the regularization parameter in BISM calculation is chosen with the chi2 criterion (see SI-1 and Fig. S6).
  * Note [2] In the experiment, the zero divergence stress tensor $\Phi$ in Eq. (5) is not apparent. This could be that the boundary conditions impose an amplitude for $\Phi$ that is buried in the noise or that its contribution is filtered by the regularization step in the calculation of $hS_{tot}$, that filters low frequencies (see SI-1).
  * Pasapera  _et al._ [2015] A. M. Pasapera, S. V. Plotnikov, R. S. Fischer, L. B. Case, T. T. Egelhoff, and C. M. Waterman, Curr. Biol. 25, 175 (2015).
  * Wolfenson  _et al._ [2016] H. Wolfenson, G. Meacci, S. Liu, M. R. Stachowiak, T. Iskratsch, S. Ghassemi, P. Roca-Cusachs, B. O’Shaughnessy, J. Hone, and M. P. Sheetz, Nat. Cell Biol. 18, 33 (2016).
  * Colin-York  _et al._ [2017] H. Colin-York, C. Eggeling, and M. Fritzsche, Nat. Protoc. 12, 783 (2017).
  * Colin-York  _et al._ [2016] H. Colin-York, D. Shrestha, J. H. Felce, D. Waithe, E. Moeendarbary, S. J. Davis, C. Eggeling, and M. Fritzsche, Nano Lett. 16, 2633 (2016).
  * Nicolas  _et al._ [2004] A. Nicolas, B. Geiger, and S. A. Safran, Proc. Natl Acad. Sci. U.S.A. 101, 12520 (2004).
  * Nicolas  _et al._ [2008] A. Nicolas, A. Besser, and S. A. Safran, Biophys. J. 95, 527 (2008).