I'm mathematically modeling yellow fever, and yet, after assembling the model, I need to find the stationary states of the model. I had to do it by hand because wolfram software could not get the necessary values.

My question is, how do I check if my stationary states are in line with my biological model?

I was told that there is a way to verify this by performing a substitution of the s values of the points within my own template code. But I have no idea how to do it.

The code used is taking forever to get back to me as this data behaves in my model

function ypsir = ypsirtry(t,y)
betah = .01;
muh = 0.01;
sigmah = 0.05;
alphah = 0.06;
deltah = 0.05; 
mui = 0.07;
mum = 0.07;
betai = 0.03;
lambdai = 0.06;
lambdah = 0.5;
k = 500;

ypsir = zeros(5,1);
ypsir(1) = sigmah*y(3) - muh*y(1) - betah*y(1)*y(5) + lambdah*(y(2)+ y(3) + 
ypsir(2) = -deltah*y(2) - y(2)*(alphah + muh) + betah*y(5)*y(1);
ypsir(3) = deltah*y(2) - sigmah*y(3);
ypsir(4) = -mui*y(4) - mum*y(4) - betai*y(2)*y(4) + lambdai*(y(4) + y(5))*(1 
- (y(4) + y(5))/k);
ypsir(5) = betai*y(2)*y(4) - mui*y(5) - mum*y(5);

Along with this calling program

to = 0;
tf =50;
yo = [99 1 0 1 3];
[t y] = ode45(@(t,y)ypsirtry(t,y),[to tf],yo); 
title('First Model of Yellow Fever')
ylabel('susceptible, infected, recovered,Mosquito Non-carrier, Carrier 

My stationary states are:

S = (1/betai)*[(muh+alphah)+(sigmah/deltah)*deltah]*(mui+mum)/betai*N*I);

I = sigmah*R/deltah;

R = ((muh/betah + lambdah)*(muh+alphah))/(sigmah+lamdah)-(muh+alphah)* 
(deltah/sigmah*deltah)-(betah/muh - lambdah)*(betai*N*I)* 

N = -k/lambdai((-lamdai+(lambdai/k))*P/-betai*I + lambdai - (2k/lambdai*P)- 

P = betai*N*I/(mui+mum); 

At where S = susceptible I = infected R = retweeted N = non-carrier mosquito P = carrier mosquito

  • $\begingroup$ I think your question is a good question, but we would need more details in order to answer them. How did you code your model? Which stationary states? (Why do you use the biomart tag?) Which data do you have to verify your models (if it is not purely theoretical)? What do you mean by "performing a schedule" ? You can edit it to include this or other information. $\endgroup$
    – llrs
    Commented Mar 16, 2018 at 15:46
  • $\begingroup$ I'm already going to make the issues, sorry for the delay! I collared biomat because I thought it would be mathematical modeling in biology. I'm new to the site and I was a little lost on which site to put. I will post the code with the respective data. But when putting it to compile it is dying a lot of time. But I will add other information. Do you know how I can put on latex around here? to make it easier to see formulas in my edits. Sorry English, it is not my mother tongue how can you perceive rs ... $\endgroup$
    – wgrt
    Commented Mar 16, 2018 at 21:09
  • $\begingroup$ Wouldn't it be sufficient to determine this by running it until the system elicits stationary behavior and then comparing the values? If you did that for a variety of values for the constants (betah and so on) then that would also be a good test. $\endgroup$
    – Devon Ryan
    Commented Mar 16, 2018 at 22:31
  • $\begingroup$ @wgrt See the tour page to learn about the site and the help center page to know how it works. How much is a lot of time: 3 days? Do you know if it can be optimized? (I don't know enough matlab to propose it myself). Also you left some questions from my previous comment unanswered $\endgroup$
    – llrs
    Commented Mar 17, 2018 at 10:25
  • $\begingroup$ @Llopis You're not going to optimize an RK45 method that's implemented in C much. If this really is slow to compute then presumably there are some very drastic changes in values (RK45 shrinks its steps in such cases). $\endgroup$
    – Devon Ryan
    Commented Mar 17, 2018 at 18:25

1 Answer 1


After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.

But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct.

Regarding the model is taking forever, I took the logistics part and ran well!


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.