6

I was reading this paper:

Global stability for an HIV/AIDS epidemic model with different latent stages and treatment

Everything is understood apart from on page 7 of the pdf (page 1486 in the document). How does the author algebraically go from the second line to the last line for the equation of $\dot V$?

I don't understand how he "generates" more terms in the last line compared to the one preceding it.

EDIT: For those who cannot access the paper

The system: \begin{align*} \dot S &= \Lambda - (\beta_1 S I_2 +\beta_2 S J )-\mu S \\ \dot I_1 &= p\beta_1 S I_2 +q\beta_2 S J +\xi_1 J -b_1 I_1\\ \dot I_2 &= (1-P)\beta_1 S I_2 +(1-q)\beta_2 S J +\epsilon I_1 +\xi_2 J -b_2 I_2\\ \dot J &= p_1 I_2 -b_3 J\\ \dot A &= p_2 J - b_4 A \end{align*} Equilibrium point: \begin{align*} S^* &= \frac{\Lambda}{\mu \mathcal{R}_0}\\ I_1^* &= \frac{1}{b_1}\left[ p \beta_1 \frac{\Lambda b_3 }{\left(\beta_1 b_3 +\beta_2 p_1 \right)J^* +\mu p_1}\right. +\left. q \beta_2 \frac{\Lambda p_1}{\left(\beta_1 b_3 +\beta_2 p_1 \right)J^* +\mu p_1}+\xi_1\right]J^*\\ I_2^* &= \frac{b_3}{p_1}J^*\\ J^*&= \frac{\mu p_1 }{\beta_1 b_3 +\beta_2 p_1}\left(\mathcal{R}_0-1\right)\\ A^*&=\frac{p_2}{b_4}J^* \end{align*}

Theorem: If $p=q$ and $\mathcal{R}_0 >1$ then the above equilibrium point is globally stable.

Proof:

Define Lyapunov function: $$V = S-S^* \ln S+ B( I_1-{I_1}^* \ln I_1) +C(I_2-{I_2}^* \ln I_2) + D( J -J^* \ln J)$$

Derivative is given by:

$$\dot V =\left(1-\frac{S^*}{S}\right) \dot S + B\left(1-\frac{{I_1}^*}{I_1} \right)\dot I_1 + C\left(1-\frac{{I_2}^*}{I_2}\right)\dot I_2+D\left(1-\frac{J^*}{J} \right)\dot J$$

The author then does substitutions i.e replaces $\Lambda$, $b_1$, $b_2$, $b_3$ by making the original system equal to $0$. He finds the constants $B$, $C$ and $D$ by killing the variable co-efficients, giving:
\begin{align*} B &= \frac{\epsilon}{\epsilon p +b_1(1-p)}\\ C&= \frac{b_1}{\epsilon p +b_1(1-p)}\\ D &= \frac{b1 b_2}{p_1[\epsilon p +b_1(1-p)]} -\frac{\beta_1 S^*}{p_1} \end{align*}

Next he does another substitution $ x=\dfrac{S}{S^*}$, $y=\dfrac{I_1}{I_1^*}$, $z=\dfrac{I_2}{I_2^*}$ and $u=\dfrac{J}{J^*}$ to which he arrives at: \begin{align*} \dot V &= -\mu S^* \frac{\left(1-x\right)^2}{x}+ \left[ \beta_1 S^* {I_2}^*+\beta_2 S^* J^*+B p\beta_1 S^* {I_2}^* + B q \beta_2 S^* J^*\right.\\ &+ B \xi_1 J^*+ C(1-p)\beta_1 S^* {I_2}^*+C(1-q)\beta_2 S^* J^*\\ &+\left. C \epsilon {I_1}^*+C \xi_2 J^*+D p_1 {I_2}^*\right]-x\left[ C(1-p)\beta_1 S^* {I_2}^* \right] -\frac{xz}{y}B p\beta_1 S^* {I_2}^*\\ & -\frac{xu}{y}B q \beta_2 S^* J^* -\frac{u}{y}B \xi_1 J^* - \frac{xu}{z}C(1-q)\beta_2 S^* J^*\\ &- \frac{y}{z}C \epsilon {I_1}^* - \frac{u}{z}C \xi_2 J^*-\frac{z}{u}D p_1 {I_2}^*\\ &- \frac{1}{x}\left[\beta_1 S^* {I_2}^*+ \beta_2 S^* J^*\right]\\ &= -\mu S^* \frac{\left(1-x\right)^2}{x} + \frac{b_1}{\epsilon p +b_1(1-p)}(1-p)\beta_1 S^* I_2^* \left(2-x-\frac{1}{x}\right)\\ &+\frac{b_1}{\epsilon p +b_1(1-p)}\xi_2 J^*\left(2-\frac{u}{z}-\frac{z}{u}\right)\\ &+\frac{b_1}{\epsilon p +b_1(1-p)}(1-q)\beta_2 S^* J^* \left(3-\frac{1}{x}-\frac{xu}{z}-\frac{z}{u}\right)\\ &+\frac{\epsilon}{\epsilon p +b_1(1-p)}p \beta_1 S^* I_2^*\left(3-\frac{1}{x}-\frac{xz}{y}-\frac{y}{z}\right)\\ &+\frac{\epsilon}{\epsilon p +b_1(1-p)} q\beta_2 S^* J^*\left(4-\frac{1}{x}-\frac{y}{z}-\frac{z}{u}-\frac{xu}{y} \right)\\ &+ \frac{\epsilon}{\epsilon p +b_1(1-p)}\xi_1 J^*\left(3-\frac{y}{z}-\frac{z}{u}-\frac{u}{y} \right) \end{align*}

I don't understand how he "generates" more terms in the last line compared to the one preceding it.

  • 2
    Could you add the relevant data to the question itself? Namely, one can't guarantee that the link doesn't break in the future. – Kwin van der Veen Dec 06 '21 at 16:52
  • @KwinvanderVeen I have added some extra data, however it is really tedious typing it all out. –  Dec 07 '21 at 13:15
  • 1
    The person that might be answering this question most likely also needs to type it out, which might also hold them back from actually writing an answer. So, if they can copy it from you/the question then it might also increase the odds of your question being answered. – Kwin van der Veen Dec 07 '21 at 13:27
  • 3
    @KwinvanderVeen Okay, I will add everything. –  Dec 07 '21 at 13:30
  • Could you also add what the second and last line for $\dot{V}$ refer too? – Kwin van der Veen Dec 07 '21 at 13:30
  • 1
    @KwinvanderVeen Its just manipulation of the derivative of the Lyapunov function –  Dec 07 '21 at 13:34
  • @Math: Transcription error: The $\xi_2 J^\star$ term on your right-hand side has an extra factor of $(1-p)$. ... If I haven't made an error, then throwing everything into Mathematica and unraveling every substitutions I could find fails to verify the equality. (There's certainly more going on than simple regrouping; only one side has $I_1^$ terms.) Assuming $p=q$ and $R_0$ is the expression under $(3.13)$ gets the difference of the expressions into a relatively* compact factored form, but it's still not obviously zero. I could be overlooking a key relation in the paper. ... What a mess! – Blue Dec 08 '21 at 17:46
  • @Math: Another transcription error: In your "equilibrium point" expression for $I_1^$, the "$,J_2^,$"s should be simply "$,J^,$"s. With that fix, Mathematica* tells me that the difference in the left- and right-hand sides of the final equality reduces further to $$\epsilon p_1\xi_1;(;\Lambda (b_3\beta_1+\beta_2p_1)H-\mu K;);(;I_1 b_1b_3H-I_2pK;);\stackrel{?}{=};0$$ where $H:=\epsilon p+b_1(1-p)$ and $K:=b_1b_2b_3-\epsilon p_1\xi_1-b_1p_1\xi_2$. Why any factor vanishes remains unclear to me. (Of course, there's still a possibility of my own error here.) – Blue Dec 09 '21 at 09:35
  • @Blue I believe there shouldn’t be any variables($S, I_1, I_2, J, A$) in the equation for $\dot V$. So your last equality seems off? Or did you mean $I_1^$ and $I_2^$? – Leo Dec 09 '21 at 10:10
  • @Leo: Well, the whole point of this question is that last equality shouldn't have any variables at all; it "should be" zero! :) ... The good news is that, after yet another pass at my Mathematica input, I did indeed get to zero, but it took back-substituting $B$, $C$, $D$, and $S^$, $I_1^$, $I_2^$, $J^$, and $R_0$, and $p=q$. (And fixing my own typo.) So, OP's concern that the target equality doesn't obviously follow algebraically appears justified. The bad news is that "Mathematica says so" isn't a satisfactory response. :) – Blue Dec 09 '21 at 10:57
  • 1
    @Blue Can you post your Mathematica code? When I tried, it wasn't equalling.. –  Dec 09 '21 at 11:31

1 Answers1

1

Posting Mathematica code, as requested ...

(Some notational conventions: I write bb for $\beta$, ee for $\epsilon$, xx for $\xi$, LL for $\Lambda$, mm for $\mu$; a trailing s corresponds to a superscript $*$, as in Ss for $S^*$; and I tend to double my uppercase variables ---BB, CC, DD for $B$, $C$, $D$--- to avoid conflicting with Mathematica's predefined C and D operators. Also, I don't bother with index subscripts; eg, bb1 is $\beta_1$ and I2s is $I_2^*$. )

The first pass at calculating the difference of the left- and right-hand sides of the target equality. I did some preliminary grouping, ignored the common first term on each side, and also temporarily re-introduced $B$ and $C$ for the right-hand side to reduce clutter):

Factor[(
  Ss (bb1 I2s + bb2 Js) (1 - 1/x)
+ BB ( p bb1 I2s Ss (1 - x z/y)
     + q bb2 Js Ss (1 - x u/y)
     + xx1 Js (1 - u/y) )
+ CC ( (1 - p) bb1 I2s Ss (1 - x) 
     + (1 - q) bb2 Js Ss (1 - x u/z) 
     + ee I1s (1 - y/z)
     + xx2 Js (1 - u/z) )
+ DD p1 I2s (1 - z/u) 
) - (
  BB ( p bb1 I2s Ss (3 - 1/x - x z/y - y/z)
     + q bb2 Js Ss (4 - 1/x - y/z - z/u - x u/y)
     + xx1 Js (3 - y/z - z/u - u/y) )
+ CC ( (1 - p) bb1 I2s Ss (2 - x - 1/x)
     + xx2 Js (2 - u/z - z/u)
     + (1 - q) bb2 Js Ss (3 - 1/x - x u/z - z/u) )
)]

Terms from the $B$ and $C$ group cancel, but there's still a bit of a mess left.

Second pass, substituting the $B$, $C$, $D$ definitions:

Factor[% /. {
         DD -> b1 b2/p1/(ee p + b1 (1 - p)) - bb1 Ss/p1,
         BB -> ee/(ee p + b1 (1 - p)),
         CC -> b1/(ee p + b1 (1 - p)) }] 

The mess expands.

Third pass, substituting the "equilibrium point" expressions, and $q\to p$:

Factor[% /. {
  Ss  -> LL/mm/R0,
  I1s -> Js/b1 ((p bb1 LL b3 + q bb2 LL p1)/((bb1 b3 + bb2 p1) Js + mm p1) + xx1),
  I2s -> b3/p1 Js } /. {
  Js  -> p1 mm/(bb1 b3 + bb2 p1) (R0 - 1),
  q   -> p}]

This gives the following result. (I'm omitting the denominator, and factors R0-1 (which the article assumes is non-zero in the given context) and u-z.)

  LL ( ee p + b1 (1-p) ) (b3 bb1 + bb2 p1)
- R0 mm ( b1 b2 b3 - p1 (ee xx1 + b1 xx2) )

The article (after equation (3.13)) calculates that, for the situation under consideration, $$R_0 = \frac{(\epsilon p + b_1 (1-p))(b_3 \beta_1 + \beta_2 p_1)\frac{\Lambda}{\mu}}{b_1b_2b_3-p_1(\epsilon \xi_1+b_1\xi_2)}$$ which is precisely what is needed to make the above expression vanish. $\square$

There may be a possibility of getting to zero without explicitly unravelling every substitution; exploring that possibility is left as an exercise to the reader.

Blue
  • 83,939
  • This is nice, but now the question arises; if you didn't have the last expression, how would you come to that without knowing it? –  Dec 09 '21 at 14:01
  • Also, why didn't this work when I put everything into Mathematica? –  Dec 09 '21 at 14:15
  • @Math: "if you didn't have the last expression, how would you come to that without knowing it?" Presumably it's explained in the article, but deciphering that stuff is beyond the scope of my interest in this problem. ;) However, what the above shows is that you can't verify the target equality without that value of $R_0$; just back-substituting $B$, $C$, $D$, $S^$, $I_1^$, $I_2^$, $J^$, and $p=q$ isn't enough. – Blue Dec 09 '21 at 14:16
  • @Math: "Also, why didn't this work when I put everything into Mathematica?" If it's any consolation, it took me a few tries to just to eliminate my own typos. :) – Blue Dec 09 '21 at 14:17
  • If I attach a slightly different model as an addendum, do you think you find the last equality? –  Dec 09 '21 at 14:30
  • @Math: "If I attach a slightly different model as an addendum, do you think you find the last equality?" If I can't, someone else might. That said, if there may be considerable work involved, it might be better to post the derivation as a separate question (and link back to this one for context). Greatly broadening the scope of a question after you have received an answer is inappropriate. – Blue Dec 09 '21 at 14:36
  • 2
    No problem, I will post a separate question and refer it to you since you have the context. –  Dec 09 '21 at 14:39
  • 2
    Please check this question: https://math.stackexchange.com/questions/4328385/how-does-this-expression-follow-algebraically-from-the-last-one-continued, I will put a bounty on it whenever I can. –  Dec 09 '21 at 14:58