2
$\begingroup$

Install the following package:

PacletInstall["EcoEvo", 
 "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]

<< EcoEvo`

The code:

SetModel[{Pop[pop] -> {
    Component[
      s] :> {Equation :> \[Nu] - \[Beta]1 i2 s  - \[Beta]2 j s - \
\[Mu] s},
    Component[
      i1] -> {Equation :>  
       p \[Beta]1 i2 s + q \[Beta]2 j s - b1 i1 + \[Xi]1 j }, 
    Component[
      i2] -> {Equation :>  (1 - p) \[Beta]1 i2 s + (1 - 
           q) \[Beta]2 j s + \[Epsilon] i1 - b2 i2 + \[Xi]2 j }, 
    Component[j] -> {Equation :> p1 i2 - b3 j },
    Component[a] -> {Equation :> p2 j  - b4 a}}, 
  Parameters :> {\[Nu] > 0, \[Beta]1 > 0, \[Beta]2 > 0, \[Mu] > 0, 
    p > 0, q > 0, b1 > 0, b2 > 0, b3 > 0, b4 > 0, 
    p1 > 0, \[Epsilon] > 0}}]
eq = SolveEcoEq[]

where

 r0 = ((\[Beta]1 b3 (\[Epsilon] p + 
          b1 (1 - p)) + \[Beta]2 p1 (\[Epsilon] q + 
          b1 (1 - q))) (\[Nu]/\[Mu]))/(b1 b2 b3 - (\[Epsilon] \[Xi]1 +
         b1 \[Xi]2) p1);

We obtain our equilibrium points:

{{s -> \[Nu]/\[Mu], i1 -> 0, i2 -> 0, j -> 0, 
  a -> 0}, {s -> (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
    b1 p1 \[Xi]2)/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
    b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
    p1 q \[Beta]2 \[Epsilon]), 
  i1 -> (-b2 b3^2 p \[Beta]1 \[Nu] - b2 b3 p1 q \[Beta]2 \[Nu] - 
      b3 p1 \[Beta]1 \[Nu] \[Xi]1 + b3 p p1 \[Beta]1 \[Nu] \[Xi]1 - 
      p1^2 \[Beta]2 \[Nu] \[Xi]1 + p1^2 q \[Beta]2 \[Nu] \[Xi]1 + 
      b3 p p1 \[Beta]1 \[Nu] \[Xi]2 + p1^2 q \[Beta]2 \[Nu] \[Xi]2 + (
      b2 b3^2 p \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b2 b3 p1 q \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b3 p1 \[Beta]1 \[Mu] \[Xi]1 (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b3 p p1 \[Beta]1 \[Mu] \[Xi]1 (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - p1 q \[Beta]2 \[Epsilon]) + (
      p1^2 \[Beta]2 \[Mu] \[Xi]1 (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      p1^2 q \[Beta]2 \[Mu] \[Xi]1 (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 +
          b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b3 p p1 \[Beta]1 \[Mu] \[Xi]2 (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - p1 q \[Beta]2 \[Epsilon]) - (
      p1^2 q \[Beta]2 \[Mu] \[Xi]2 (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 +
          b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]))/(-b1 b2 b3^2 \[Beta]1 - 
      b1 b2 b3 p1 \[Beta]2 + b3 p1 \[Beta]1 \[Epsilon] \[Xi]1 + 
      p1^2 \[Beta]2 \[Epsilon] \[Xi]1 + b1 b3 p1 \[Beta]1 \[Xi]2 + 
      b1 p1^2 \[Beta]2 \[Xi]2), 
  i2 -> (-b1 b3^2 \[Beta]1 \[Nu] + b1 b3^2 p \[Beta]1 \[Nu] - 
      b1 b3 p1 \[Beta]2 \[Nu] + b1 b3 p1 q \[Beta]2 \[Nu] - 
      b3^2 p \[Beta]1 \[Epsilon] \[Nu] - 
      b3 p1 q \[Beta]2 \[Epsilon] \[Nu] + (
      b1 b3^2 \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b1 b3^2 p \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b1 b3 p1 \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b1 b3 p1 q \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b3^2 p \[Beta]1 \[Epsilon] \[Mu] (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - p1 q \[Beta]2 \[Epsilon]) + (
      b3 p1 q \[Beta]2 \[Epsilon] \[Mu] (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]))/(-b1 b2 b3^2 \[Beta]1 - 
      b1 b2 b3 p1 \[Beta]2 + b3 p1 \[Beta]1 \[Epsilon] \[Xi]1 + 
      p1^2 \[Beta]2 \[Epsilon] \[Xi]1 + b1 b3 p1 \[Beta]1 \[Xi]2 + 
      b1 p1^2 \[Beta]2 \[Xi]2), 
  j -> (-b1 b3 p1 \[Beta]1 \[Nu] + b1 b3 p p1 \[Beta]1 \[Nu] - 
      b1 p1^2 \[Beta]2 \[Nu] + b1 p1^2 q \[Beta]2 \[Nu] - 
      b3 p p1 \[Beta]1 \[Epsilon] \[Nu] - 
      p1^2 q \[Beta]2 \[Epsilon] \[Nu] + (
      b1 b3 p1 \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b1 b3 p p1 \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b1 p1^2 \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b1 p1^2 q \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b3 p p1 \[Beta]1 \[Epsilon] \[Mu] (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - p1 q \[Beta]2 \[Epsilon]) + (
      p1^2 q \[Beta]2 \[Epsilon] \[Mu] (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]))/(-b1 b2 b3^2 \[Beta]1 - 
      b1 b2 b3 p1 \[Beta]2 + b3 p1 \[Beta]1 \[Epsilon] \[Xi]1 + 
      p1^2 \[Beta]2 \[Epsilon] \[Xi]1 + b1 b3 p1 \[Beta]1 \[Xi]2 + 
      b1 p1^2 \[Beta]2 \[Xi]2), 
  a -> (-b1 b3 p1 p2 \[Beta]1 \[Nu] + b1 b3 p p1 p2 \[Beta]1 \[Nu] - 
      b1 p1^2 p2 \[Beta]2 \[Nu] + b1 p1^2 p2 q \[Beta]2 \[Nu] - 
      b3 p p1 p2 \[Beta]1 \[Epsilon] \[Nu] - 
      p1^2 p2 q \[Beta]2 \[Epsilon] \[Nu] + (
      b1 b3 p1 p2 \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b1 b3 p p1 p2 \[Beta]1 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 +
          b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b1 p1^2 p2 \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) - (
      b1 p1^2 p2 q \[Beta]2 \[Mu] (-b1 b2 b3 + p1 \[Epsilon] \[Xi]1 + 
         b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + b1 b3 p \[Beta]1 - 
       b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]) + (
      b3 p p1 p2 \[Beta]1 \[Epsilon] \[Mu] (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - p1 q \[Beta]2 \[Epsilon]) + (
      p1^2 p2 q \[Beta]2 \[Epsilon] \[Mu] (-b1 b2 b3 + 
         p1 \[Epsilon] \[Xi]1 + b1 p1 \[Xi]2))/(-b1 b3 \[Beta]1 + 
       b1 b3 p \[Beta]1 - b1 p1 \[Beta]2 + b1 p1 q \[Beta]2 - 
       b3 p \[Beta]1 \[Epsilon] - 
       p1 q \[Beta]2 \[Epsilon]))/(-b1 b2 b3^2 b4 \[Beta]1 - 
      b1 b2 b3 b4 p1 \[Beta]2 + b3 b4 p1 \[Beta]1 \[Epsilon] \[Xi]1 + 
      b4 p1^2 \[Beta]2 \[Epsilon] \[Xi]1 + 
      b1 b3 b4 p1 \[Beta]1 \[Xi]2 + b1 b4 p1^2 \[Beta]2 \[Xi]2)}}

These are long and tedious. Is there a way to simplify them using each other and r0?

For example; why does i2 = $\frac{b_3}{p_1} j$?

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ Is this a question about Mathematica? or the Wolfram Language? If you include the points (or representative data) and the Mma code that you have tried out so far for the simplification that you have alluded to, then this question will likely be considered valid and on-topic given that you can also describe where you are having difficulty. As posed, this question is asking the potential respondents to read a scientific paper. I would request you to edit your post further. $\endgroup$
    – Syed
    Commented Nov 15, 2021 at 12:52
  • $\begingroup$ @Syed I edited the question. $\endgroup$
    – Math
    Commented Nov 15, 2021 at 12:58
  • 1
    $\begingroup$ The situation has almost not changed. Now the respondents are being asked to install and understand the workings of a package related to the scientific paper. Maybe someone who has prior experience with this package would take up this question. You will have to wait and see. Why not contact the author of the paper directly? Best of luck. $\endgroup$
    – Syed
    Commented Nov 15, 2021 at 13:51
  • $\begingroup$ It isn't a difficult task to install the package to be honest. You just copy, paste and shift+ enter. All I ask is for example why is i2, obtained from our equilibrium point equal to $\frac{b_3}{p_1} j$? What simplification can we do to achieve this result? $\endgroup$
    – Math
    Commented Nov 15, 2021 at 13:54
  • $\begingroup$ @rhermans You don't require the link to answer the question since I edited the question. Let me further edit it. Please have a look now. I linked the paper as a reference in case one asked "how did you get this simplified expression?" $\endgroup$
    – Math
    Commented Nov 15, 2021 at 14:25

2 Answers 2

Reset to default
5
$\begingroup$

eqPts is the equilibrium points in your question. Start by using Simplify

LeafCount /@ {eqPts, eqPts2 = eqPts // Simplify}

(* {1958, 448} *)

The variables are

vars = Variables[Level[eqPts2, {-1}]]

(* {a, b1, b2, b3, b4, i1, i2, j, p, p1, p2, q, s, β1, β2, \
ϵ, μ, ν, ξ1, ξ2} *)

Convert the replacement rules in eqPts2 into the corresponding equations

eqns = eqPts2 /. Rule :> Equal;

With five equations you can Solve for one variable while eliminating four variables.

sol = Solve[#, i2, {s, i1, a, b1}] & /@ eqns

(* {{}, {{i2 -> (b3 j)/p1}}} *)

The expected result,

sol[[2, 1, 1]] /. Rule :> Equal

(* i2 == (b3 j)/p1 *)

Verifying this relation at both equilibrium points,

% /. eqPts2 // Simplify

(* {True, True} *)

EDIT: A more general approach will be much slower

vars = eqns[[1, All, 1]]

(* {s, i1, i2, j, a} *)

params = Complement[Variables[Level[eqns, {-1}]], vars]

(* {b1, b2, b3, b4, p, p1, p2, q, β1, β2, ϵ, μ, ν, ξ1, ξ2} *)

solve[solveVar_Symbol, var_Symbol] :=
 SortBy[
   Union[
    Simplify[
     Solve[eqns[[2]], solveVar, #] & /@
       (Append[Complement[vars,
            {solveVar, var}], #] & /@
         params) /. {} :> Nothing]],
   LeafCount][[1, 1]]

To Solve for i2 in terms of j

solve[i2, j]

(* {i2 -> (b3 j)/p1} *)

To Solve for i2 in terms of j

solve[i1, j] // Apart

(* {i1 -> (j (b3 p β1 + p1 q β2) ν)/(
   b1 (b3 j β1 + j p1 β2 + p1 μ)) + (j ξ1)/b1} *)

You could then generalize this to solve for a variable in terms of each of the other variables and select the simplest form; however, this would be very slow.

Improve this answer
$\endgroup$
10
  • $\begingroup$ Thank you very much. Can I edit the extra parts(the other components) to see if you can simplify them too? $\endgroup$
    – Math
    Commented Nov 15, 2021 at 16:50
  • $\begingroup$ You now have an approach, try to apply it to your other relations. If you encounter a problem that you cannot solve with reasonable effort, post a new question showing what you have tried and the problems that you are having. $\endgroup$
    – Bob Hanlon
    Commented Nov 15, 2021 at 16:55
  • $\begingroup$ I will try, thank you. $\endgroup$
    – Math
    Commented Nov 15, 2021 at 16:56
  • $\begingroup$ I forgot to ask; will this work without knowing the "simplified expression" beforehand? So a better example would've been the one I commented on Chris' answer. $\endgroup$
    – Math
    Commented Nov 16, 2021 at 11:24
  • $\begingroup$ I tried your edited version, it doesn't seem to work.. at least on my computer. can you post the full code? $\endgroup$
    – Math
    Commented Nov 18, 2021 at 13:44
3
$\begingroup$

$i_2 = \frac{b_3}{p_1} j$ follows directly from setting ${dj \over dt}=0$. You could do that manually, or in EcoEvo, give SolveEcoEq a list of variables to solve for and use QSS -> True (otherwise non-specified variables will be assumed equal to zero).

SolveEcoEq[{j}, QSS -> True]
(* {{j -> (i2 p1)/b3}} *)

You might also try:

SolveEcoEq[{a}, QSS -> True]
(* {{a -> (j p2)/b4}} *)

The others don't seem particularly useful.

Sorry, I don't know of any easy way to rearrange in terms of r0.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I completely overlooked the equation $j$ and this would give the result directly as stated. Now however, a better example would be why does $$ i_1 = \frac{1}{b_1}\left[p \beta_1 \frac{\nu b_3}{(\beta_1 b_3+\beta_2 p_1)J^*+\mu p_1}+q \beta_2 \frac{\nu p_1}{(\beta_1 b_3+\beta_2 p_1)J^*+\mu p_1} +\xi_1\right]J^*$$ $\endgroup$
    – Math
    Commented Nov 16, 2021 at 11:15
  • $\begingroup$ Mathematica's good at solving equations, but not necessarily putting the solution in the form that makes intuitive sense. That might still require human intervention! $\endgroup$
    – Chris K
    Commented Nov 16, 2021 at 13:54
  • $\begingroup$ Essentially we are selecting one component and describing the rest using this component, I think Mathematica should be able to do it! if you have a look at this question: mathematica.stackexchange.com/questions/258344/… we see most of components are in terms of the $J^*$ component. $\endgroup$
    – Math
    Commented Nov 16, 2021 at 14:06
  • $\begingroup$ I'll be curious to see what other people come up with. Anyhow, sometimes you gain more intuitive understanding about a model by pushing through the algebra by hand, rather than relying on Mathematica. I recommend using the two together. $\endgroup$
    – Chris K
    Commented Nov 16, 2021 at 14:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.