(* Dr. Jessica Clark University of Edinburgh & University of Glasgow*)

(*MS title - Pathogen Dynamics Across the Diversity of Ageing*)

(*Analytical analysis of the effects of individual senescence on population dynamics*)

(*uninfected population dynamics*)

ClearAll["Global`*"]

n = Yyu+You + Oyu + Oou

dYyu= (r*(1+x)*Yyu+r*(1+x)*You)*(1-(n/k))-((d*(1-y))+m)*Yyu

dYou=(r*(1-x)*Oyu+r*(1-x)*Oou)*(1-(n/k)) - ((d*(1-y))+m)*You

dOyu=m*Yyu-((d*(1+y))*Oyu)

dOou= m*You-((d*(1+y))*Oou)

Sol = FullSimplify[Solve[{dYyu==0&&dYou==0&&dOyu==0&&dOou ==0},{Yyu,You, Oyu,Oou}]]

Set @@@ Sol[[2]]

total = Yyu + You + Oyu + Oou

FullSimplify[total]

FullSimplify[D[total,x]]

x=0;y=0;FullSimplify[-((d k (1+y) (d+m-d y) (m-d (1+y)))/(r (((m-m x+d (1+x) (1+y)))^2)))]

FullSimplify[Reduce[((d k (d-m))/((d+m) r))>0&&d>0&&m>0&&k>0&&y≥-1&&y≤1&&x≥-1&&x≤1&&r>0]]

ClearAll[x,y];FullSimplify[D[total,y]]

x=0;y=0;FullSimplify[((d k ((m^2) (-1+x)-2 d m (-1+x) y+(d^2) (1+x) (((1+y))^2)))/(r (((m-m x+d (1+x) (1+y)))^2)))]

FullSimplify[Reduce[((d k (d-m))/((d+m) r))>0&&d>0&&m>0&&k>0&&y≥-1&&y≤1&&x≥-1&&x≤1&&r>0,m]]

(*------------------------------------------------------------------*)

(*------------------------------------------------------------------*)

(*------------------------------------------------------------------*)

(*Next generation method*)

ClearAll["Global`*"]

(*making the F (new infection terms), V- (infection class losses) and V+ (infection class gains) vectors*)

Ni =Yyi+Yoi+Oyi+Ooi

Byy = B*(1-A)*(1-M)

Byo=B*(1-A)*(1+M)

Boy=B*(1+A)(1-M)

Boo=B*(1+A)(1+M)

F={(Byy*Yyu*Ni),(Byo*You*Ni),(Boy*Oyu*Ni),(Boo*Oou*Ni)}
MatrixForm[F]

Vm={(Yyi*(m+v+c+(d*(1-y)))),(Yoi*(m+v+c+(d*(1-y)))), (Oyi*(v+c+(d*(1+y)))),(Ooi*(v+c+(d*(1+y))))}

MatrixForm[Vm]

Vp={0,0,m*Yyi, m*Yoi}

MatrixForm[Vp]

(*subtract the minuses from the pluses*)

V=Vm-Vp

MatrixForm[V]

(*jacobian matrices of the f and v matrices *)
(*the v matrix is made weirdly and called insandouts because mathematica didn't like it being called v or anything with v in it and made it hard to make*)

v11=D[V[[1]],Yyi]

v12=D[V[[1]],Yoi]

v13=D[V[[1]],Oyi]

v14=D[V[[1]],Ooi]

v21=D[V[[2]],Yyi]

v22=D[V[[2]],Yoi]

v23=D[V[[2]],Oyi]

v24=D[V[[2]],Ooi]

v31=D[V[[3]],Yyi]

v32=D[V[[3]],Yoi]

v33=D[V[[3]],Oyi]

v34=D[V[[3]],Ooi]

v41=D[V[[4]],Yyi]

v42=D[V[[4]],Yoi]

v43=D[V[[4]], Oyi]

v44=D[V[[4]],Ooi]

insandouts={{v11,v12,v13,v14},{v21,v22,v23,v24},{v31,v32,v33,v34},{v41,v42,v43,v44}}

f={{D[F[[1]],Yyi],D[F[[1]],Yoi],D[F[[1]],Oyi],D[F[[1]],Ooi]},{D[F[[2]],Yyi],D[F[[2]],Yoi],D[F[[2]],Oyi],D[F[[2]],Ooi]},{D[F[[3]],Yyi],D[F[[3]],Yoi],D[F[[3]],Oyi],D[F[[3]],Ooi]}, {D[F[[4]],Yyi],D[F[[4]],Yoi],D[F[[4]],Oyi],D[F[[4]],Ooi]}}

(*evaluate at the disease free equilibrium*)

n = Yyu+You+Oyu+Oou

dYyu= (r*(1+x)*Yyu+r*(1+x)*You)*(1-(n/k))-((d*(1-y))+m)*Yyu

dYou=(r*(1-x)*Oyu+r*(1-x)*Oou)*(1-(n/k)) - ((d*(1-y))+m)*You

dOyu=m*Yyu-((d*(1+y))*Oyu)

dOou= m*You-((d*(1+y))*Oou)

Sol = FullSimplify[Solve[{dYyu==0&&dYou==0&&dOyu==0&&dOou ==0},{Yyu,You, Oyu,Oou}]]

Set @@@ Sol[[2]]

f

insandouts

(*create the next generation matrix*)

output=f.Inverse[insandouts]

(*the dominant eigenvalue is the equivalent of R0*)

FullSimplify[Eigenvalues[output]]

R0=(B k (((-1+A) d-(1+A) m) (c+d+m+v)+d ((-1+A) c+2 (-1+A) d+2 A m+(-1+A) v) y+(-1+A) (d^2) (y^2)) (m (1+M) (-1+x)+d (-1+M) (1+x) (1+y)) (-m r (-1+x)-d (m-r (1+x)) (1+y)+(d^2) (-1+(y^2))))/(r (c+d+m+v-d y) (d+m+d y) (c+d+v+d y) (((m-m x+d (1+x) (1+y)))^2))

(*R0*)

ClearAll[A,M,B,x,y,d,r,m,k]

Dr0x = D[R0,x]

y = 0; x=0; A= 0; M =0

(*DR0X*)

FullSimplify[Dr0x]

Simplify[Reduce[Dr0x >0 &&B>0&&d>0&&k>0&&r>0&& m >0&&v>0&&c>0]]

(* Dr0y*)

ClearAll[y,A,M,x];Dr0y = D[R0,y]

y = 0;x=0;A=0;M=0;FullSimplify[Dr0y]

FullSimplify[Reduce[Dr0y >0 &&B>0&&d>0&&k>0&&r>0&& m >0&&v>0&&c>0,d]]

(*Dr0y*)

ClearAll[y,A,M,x];Dr0A = D[R0,A]

y = 0;x=0;A=0;M=0;FullSimplify[Dr0A]

FullSimplify[Reduce[Dr0A >0 &&B>0&&d>0&&k>0&&r>0&& m >0&&v>0&&c>0&&r>d,r]]