Message Boards Message Boards


[WSG21] Daily Study Group: Building and Applying Epidemic Models

Posted 3 months ago
52 Replies
14 Total Likes

A new study group for Building and Applying Epidemic Models with the Wolfram Language begins Monday, May 10, 2021!

Making progress in an online course can be daunting when you have to study all alone. Join a cohort of fellow Wolfram Language users for a two-week study group in which you will start in week one with the basics of implementing compartment-based epidemiological models in the Wolfram Language. In the second week, you will cover multi-group models–taking into account vital demographics (i.e. birth, death) and age groups–in addition to introducing control measures and ending with stochastic models.

Sign up here:

52 Replies
Posted 3 months ago


The c19Data = ResourceFunction["JHUCOVID19VaccineData"][]

will not load on my machine running Mathematica Desktop. It runs forever then gives a kernal dead error .

When I run it in the cloud , it runs but gives a number of warnings

enter image description here

Posted 2 months ago

OK , it seems the menu

Evaluation \Dynamic Updating Enabled

must be enabled (have a tick )

then it loads

Related to the question asked, I would like to know the following: From c19Data = ResourceFunction["JHUCOVID19VaccineData"][] how do I determine what the keys are, and how do I extract a set of values, for example the time series corresponding to TotalVaccineDosesAllocated?

Posted 2 months ago

I used

Keys[c19Data ][[1]]

c19Data [[All,{"AdministrativeDivision","TotalVaccineDosesAllocated"}]]

Thank you. It worked well. If I may trouble you again: I did following your instructions for the ResourceFunction["COVID19EpidemicData"] found in WSG14_Session1.nb but got an unexpected result (as shown in the attachment). What did do wrong?


Posted 2 months ago

Thank you.

Sir, this might be silly and I am kinda new to Mathematica. I have been a MATLAB and Python user always. When using ParametricNDSolve to solve for solving and plotting, it gives the plot of all the variables S, E, I, and R. I want to vary the parameter using the slider but only focus on the Infection and Recovered compartment. I tried doing it but it showed an error. I am showing what I tried here in the attached notebook. Thank You so much.

Posted 2 months ago

Hi Ravi,

It is easier to pick the solutions you want to plot from the full list. They are in S, E, I, R order, I and R are the third and fourth.

With[{seirSolutions = 
   seirCompartments /. 
    ParametricNDSolve[seirEquations /. seirInitialConditions, seirCompartments, {t, 0, 100000},
    {\[Beta], \[Gamma], \[Sigma]}]}, 
      Through[seirSolutions[\[Beta], \[Gamma], \[Sigma]]][[3 ;; 4]], {t, 0, tmax}, 
     PlotRange -> {{0, tmax}, {0, ymax}}]}],
     (*paramters in model*)
     {{\[Beta], 1/2}, 0, 1}, {{\[Gamma], 1/10}, 0, 1}, {{\[Sigma], 3/4}, 0, 1}, 
    (*plot range limits*)
   {{tmax, 100, "\!\(\*SubscriptBox[\(t\), \(max\)]\)"}, 0, 200}, 
   {{ymax, 1000000, "\!\(\*SubscriptBox[\(y\), \(max\)]\)"}, 0, 1000000}, 
  SaveDefinitions -> True]]

enter image description here

Thank you so much, Rohit. This works if they are in order right. can I plot if they are not in order, say S and R.

Edit: Okay I got this from wolfram documentation on Part if anyone needs it here it is

NDSolve and ParametricNDSolve both return lists of Rules. I find it more convenient to work with those lists of Rules than just a list of their Values (right hand sides) as is done in your code. This way I work directly with the variables I need rather than some part of a list. (In general I try to avoid using Part.)

I also find it easier to solve for the variables vi rather than vi[t].

I made a couple changes to your code. In the first cell I used

seirCompartments = {\[ScriptCapitalS], \[ScriptCapitalE], \
\[ScriptCapitalI], \[ScriptCapitalR]};

and in the second cell I used

Plot[{\[ScriptCapitalS][\[Beta], \[Gamma], \[Sigma]][
     t], \[ScriptCapitalR][\[Beta], \[Gamma], \[Sigma]][t]} /. 
   seirSolutions // Evaluate, {t, 0, tmax}, 
 PlotRange -> {{0, tmax}, {0, ymax}}]
Posted 2 months ago

Where are the study notebooks, please? I have joined very late....

Register here

Then you should be able access the the materials.

Posted 2 months ago

Post registration (now done 3 times) I land here: but I see no NB links.

The series materials are provided daily in the invitation link you get to join the session. You will get them there, but for now, if you have registered here is the link where you can get notebooks. Files for the series

Hopefully, I am not doing something wrong by sharing this.

Posted 2 months ago


Posted 2 months ago

In the DynamicTransmissionModel function, if the option AgeStratification is not enabled, the ForceOfInfection is effectively the probability of one of the contacts made by the suseptibles being with one individual in an infective compartment; if the AgeStratification is enabled, the ForceOfInfection is effectively the multiplication of that probability and the contact rates. Why does this inconsistence exist?

Regarding the force of infection in age-stratified vs unstratified models, the really isn't an omission of the contract rate from the unstratified case. When one writes the transition

Transition[S, \[Beta] \[Lambda][t], I]

the interpretation of the parameter \[Beta] is up to the user, and as was mentioned in the sessions it is often a composition of several effects including probability of transmission and contact rate. If the code inserted a contact rate parameter in the force of infection for the unstratified case

\[Lambda][t] -> m I[t]

then the infection term in the equations would be

\[Beta] m S[t] I[t]

The two parameters would always appear together and a single composite parameter would suffice.

In the age stratified case, a contact matrix is used and the force of infection is

Subscript[\[Lambda], i][t] -> Sum[Subscript[M, i, j] Subscript[I, j][t], {j, nAges}]

and the infection term in the equations would be

\[Beta] Subscript[S, i][t] Sum[Subscript[M, i, j] Subscript[I, j][t], {i, nAges}]

The parameter \[Beta] would appear with different elements of the contact matrix, and a single parameter does not suffice.

It's really a matter of detail about how one wants to model a combination of effects, and of convenience.

Posted 2 months ago

Has the Stochastic Models video been released ?

Posted 2 months ago

Can you give any analytical insight into "infections/lives saved" as function of vaccine rollout rate? Clearly it is easy, if vaccination is "instant" and I have had success assuming vaccine rollout is far faster than new infections. I appreciate that "lives" is harder that "cases", with age stratification on both vaccine rollout and mortality but that is not my question.

Computing the number of lives saved, money saved, quality adjusted life years (QALYs) gained, etc. are one of the leading reasons pharmaceutical companies build these epidemiology models. In a nutshell, the model is run without vaccination as the base case, and then for several vaccination strategies (varying age of vaccination, number of doses, assumptions about efficacy, etc.). The health and economics outcomes are computed for each strategy (including no-vaccine) and compared to each other. The results are used to choose the best strategy to propose to regulators (FDA).

The strategies are typically simulated over a 10, 25, 50, or 100 year time horizon because the results are often time-dependent. Sometimes it takes several years for the vaccination program to reach its full effect in a population.

Posted 2 months ago

I was thinking of Covid19 rollouts, key decisions which need to be made on a much shorter time horizon (e.g. the UK decision to prioritorize first doses and older people). But even then I appreciate the “best” simulation (e.g. with full age stratification ..) is needed.

But my interest is the insight one can get from the simplest models which, in select cases, can even be done on the “back of an envelope”. I have attached and example Notebook, with no “progressing” vaccination.

I am extending for the actual rollout but, so far, it is not accurate rate and/or elegant enough.

Posted 2 months ago

If I want the absolute vaccination rate to be controlled do I really need a "negative sink". see attached .nb

Posted 2 months ago

It would be excellent if you (Robert N) could comment on this some time (it relates to he previous post)

Great question! Actually, you implemented a "negative source" in the model, and if you look at the model equations generated by DynamicTransmissionModel, they are just what you want. You do have to be careful that you don't consume all the susceptible individuals and S[t] < 0 during the simulation.

I admit that it's a bit awkward to do it this way. In systems biology, SBML specifies the models in a very similar way to CompartmentalModeling.wl. However, SBML also permits the user to specify a kinetic rate law for each transition which allows one to use 0-order kinetics (like you want here), first- and second-order (mass action), or catalyzed/inhibited (Michaelis-Menten or Hill). It's in my plans to add this capability to CompartmentalModeling.wl.

I'll comment later about your questions/remarks highlighted in yellow, I just need some time to go through the notebook carefully.

Posted 2 months ago

Many thanks for your responses.

It did, indeed, make sure I did not use up all my susceptible. For our discussions I removed another detail: in many countries we vaccinate irrespective of antibodies/ past infections, so needing a factor of approx. (1-R[t]).

While I was a Control Engineer decades ago, for this I am a “hobby epidemiologist. I have done loads of SIR modelling and am well familiar with the herd immunity that one can achieve via recovered cases and vaccination. E.g. effective reproduction number

Re[t] = k Ro S[t] = k Ro (1 – R[t] – V[t])) = 1 for herd immunity, where k is for NPIs or season. They seem like disease free equilibria to me.

But I am unclear if these are valid equilibria for NextGenerationMatrix and, if so, what parameters to use to extract them. When I flip from “your” vaccine model with S[t]*phi[t] to mine, with vt[t]. none are found. I note yours is inextricably linked to VitalDemographics while mine is not.

Which is why I sent my “too long” post /Notebook

Is there a way to only display a short summary so Notebooks do not swamp the Forum?

Posted 2 months ago

Hello All!

Sorry for putting this up here so late, but I noticed an error I've been getting with the first notebook of this series. I don't see a population key for AdmistrativeRegion so perhaps that is why i'm getting the error so these lines of code don't render?

Error I'm Receiving

Posted 2 months ago

In Host-Vector Models, Malaria is not a bacterium but a protozoa

Thank you, I'll make the correction.

Posted 2 months ago

Dear Dr. Robert, Greetings... I am one of the attendees of this webinar. I wonder for my attached file that I couldn't get output for the DynamicTransmissionModel and consequently for the next generation matrix function? I added vital demography natural birth, natural death, and fatality death for symptomatic (ds) and asymptomatic (da) as well. But, unfortunately, I couldn't find them on the graph? Thanks for your help and support.


It appears that you are using \[RightArrow], or <esc> -><esc>, instead of \[Rule], or ->, for the options to the functions, and that that's why the functions are not evaluating. There a great many instances of this problem in the notebook, and they have to be fixed one-by-one because using Edit > Find with ReplaceAll will also change the transitions that do use the right arrow character.

The format for the values of some of the options for DynamicTransmissionModel are a bit complex, and I'm addressing that by writing much better documentation in symbol reference pages, but they are not ready yet. The option

ForceOfInfection -> {\[Lambda], Subscript[\[ScriptCapitalI], a], 
  Subscript[\[ScriptCapitalI], s] {"Standard", Automatic}}

is missing a pair of curly braces and a comma, and should be coded this way

ForceOfInfection -> {\[Lambda], {Subscript[\[ScriptCapitalI], a], 
   Subscript[\[ScriptCapitalI], s]}, {"Standard", Automatic}}

I pushed a new version of both packages to the Cloud over the weekend, and one of the changes renamed the CompartmentalModelGraph option VertexColor to CompartmentColors to avoid confusion with the Polygon option VertexColors.

Posted 1 month ago

Really many thanks. It works. But unfortunately, once I tried to include the DynamicTransmissionModel to include Birth and death rates I faced the same issue! I tried to solve but couldn't? I avoided the subscripts to avoid coding error. THANKS here is my code

There are two errors in the input:

vitalModelData = 
 DynamicTransmissionModel[model, t, 
  ForceOfInfection -> {\[Lambda], \[ScriptCapitalA], \
\[ScriptCapitalI], {"Standard", Automatic}}, 
  BirthCompartments \[RightArrow] \[ScriptCapitalS]]

The first is that the infectious compartments need to be in a list by themselves, and the second is that there was one more \[RightArrow] in the last option. This code works:

vitalModelData = 
 DynamicTransmissionModel[model, t, 
  ForceOfInfection -> {\[Lambda], {\[ScriptCapitalA], \
\[ScriptCapitalI]}, {"Standard", Automatic}}, 
  BirthCompartments -> \[ScriptCapitalS]]

By the way, are you working in a notebook on the desktop or in the Wolfram Cloud? I'm wondering if there is something I should do to make it just as easy to work in either environment.

Posted 1 month ago

Thanks for your feedback. I am working on a notebook on the desktop, not in the Wolfram Cloud. Is it better to use the Wolfram cloud? If so, I will do that. Regards

Posted 1 month ago

Unfortunately, same error, I did to group the infectious compartments together! is there any other error?



The option for the birth compartments is still using \[RightArrow] instead of \[Rule]. If you add //Hold//FullForm to the end of the line and evaluate it you can see that \[RightArrow] is getting interpreted as Transition:

In[18]:= DynamicTransmissionModel[model, t, ForceOfInfection -> {\[Lambda], 
{\[ScriptCapitalA], \[ScriptCapitalI]}, {"Standard", Automatic}}, 
BirthCompartments \[RightArrow] \[ScriptCapitalS]] // Hold // FullForm

Out[18]//FullForm= Hold[DynamicTransmissionModel[model, t, Rule[ForceOfInfection, 
List[\[Lambda], List[\[ScriptCapitalA], \[ScriptCapitalI]], List["Standard", Automatic]]], 
Transition[BirthCompartments, \[ScriptCapitalS]]]]

Regarding the working environment, I find using the desktop much more convenient that the Cloud, specially when I use a lot of typeset expressions for input.

Posted 1 month ago

Dear Dr. Nachbar, Greetings. Thanks for your continuous help and support, for my model SVEAIRS including vital demographic and vaccination, I wonder regarding the treatment model data give an output percentage as shown? moreover, the ngm command gives me an error as well? any comments or suggestions? Thanks again and regards

Posted 2 months ago

Is there a way I can default any .nb attachments to "not preview" to avoid cluttering the Forum? Or should they either be shorter, not sent and sent elsewhere?

Yes Peter, you can attach a notebook file without preview using the button "Add a file to this post" located at the end of the editor as shown in the below image. You can also publish a notebook in cloud and provide the link. enter image description here

Posted 2 months ago

Sorry, I thought that is how I had attached my two .nb that are cluttering the forum

The thread is fast and works fine with me. Feel free to use the feature that suits you best.

The ForceOfInfection option in DynamicTransmissionModel is using \[RightArrow] instead of \[Rule] or ->. This code works:

treatmentModelData = 
  DynamicTransmissionModel[treatmentModel, t, 
   ForceOfInfection -> {{\[Lambda], {\[ScriptCapitalA], \
\[ScriptCapitalI]}, {"Standard", Automatic}}}];

In[15]:= Keys@treatmentModelData

Out[15]= {"Variables", "ComponentTransitions", \
"StoichiometricMatrix", "Rates", "Equations", "TimeSymbol", \

The message from NextGenerationMatrix means that the disease-free equilibrium that was supplied does not give 0 on the rhs of the equations when the DFE is substituted in.

Looking at the model, I see that there are two places where a space for implicit multiplication is missing:

Transition[\[ScriptCapitalS], \[Beta] \[Lambda][
    t] + \[Beta] \[Omega]\[Lambda][t], \[ScriptCapitalE]]

should have a space between omega and lambda:

Transition[\[ScriptCapitalS], \[Beta] \[Lambda][
    t] + \[Beta] \[Omega] \[Lambda][t], \[ScriptCapitalE]]


Transition[\[ScriptCapitalE], \[Gamma]\[Rho], \[ScriptCapitalI]]

should have a space between gamma and rho:

Transition[\[ScriptCapitalE], \[Gamma] \[Rho], \[ScriptCapitalI]]

A space for implicit multiplication is also missing in the numerator for V of the DFE:

V -> (\[Kappa]\[CapitalLambda])/(\[Mu] (\[Kappa] + \[Mu]))

should have a space between kappa and Lambda:

V -> (\[Kappa] \[CapitalLambda])/(\[Mu] (\[Kappa] + \[Mu]))

Now, looking more closely at the model, it's not clear why the transition from S to E has two forces of infection. Moreover, the transition from V to E is missing a force of infection. After you make the corrections above, please think about your model and modify the infection transitions accordingly.

Posted 1 month ago

Really many thanks for your valuable comments. I did all the spaces as required. I wonder regarding the point of DFE where we solve all differential equations to get RHS zero. Are there any commands to find the point of DFE in wolfram?? or it is better to solve them manually (sometimes it will be very complicated to find them in a symbolic computation) any advice to evaluate them? I did the transition from S to E to have two forced due to reduction in infection based on the state of the infected individual (Asym. or Inf.) and added as well the force of infection for Vto E as well. I attached my file for your reference. Thanks again for your help and advice.


The Automatic method used in the NextGenerationMatrix function sets the derivatives in the differential equations to 0 and then solves for the model variables. This gives all the equilibria, and then the one with the infectious compartments equal to 0 is selected.

It works quite well for simple models, but can take a very long time for larger models. That's why I implemented the "DiseaseFreeEquilibrium" option. I need to find a more efficient automatic method, and one idea I had was to set the force of infection to 0 and then solve for the equilibria in the usual way.

When the state fo the infectious individual affects the transmission rate, then one still uses a single force of infection in the model and also specifies different infectivities for the infectious compartments. So, for your model use this transition:

Transition[\[ScriptCapitalS], \[Beta] \[Lambda][t], \[ScriptCapitalE]]

and this option:

ForceOfInfection -> {{\[Lambda], {{\[ScriptCapitalA], q}, {\[ScriptCapitalI], 1}}, 
{"Standard", Automatic}}}

The resulting force of infection is:

{\[Lambda][t] -> (q \[ScriptCapitalA][t] + \[ScriptCapitalI][t])/
(V[t] + \[ScriptCapitalA][t] + \[ScriptCapitalE][t] + \[ScriptCapitalI][t] + \[ScriptCapitalR][t] + 

which I believe is what you are after.

Posted 1 month ago

Is there any way to compute the endemic equilibrium points other than this code? My Wolfram Mathematica keeps running forever without any output? I wonder if there is any other solution to do so and so I will continue my analysis. Thanks in advance


That's the correct method. You could split it into two parts, one for Solve and another for Simplify. I suspect that the simplification is taking all the time, and by splitting the work you at least get the raw solutions (there may be more than one) for endemic equilibria.

It doesn't take much to make a model complicated enough to make symbolic solutions prohibitively expensive.

You might try searching some math forums (e.g., StackExchange) for methods to solve systems of simultaneous equations.

Btw, thank you for your perseverance on this problem, and for working with my package. I appreciate the dialog and valuable feedback you are giving me.

Posted 1 month ago

Thanks a lot for your clarifications, recommendations, and prompt response. Regards

Posted 1 month ago

Dear Dr. Robert, I am doing my covid19 model related to UAE and am interested to complete it using your amazing tool and techniques. I have read all your previous models as well and have a doubt regarding your "Epidemiological Models for Influenza and COVID-19", how we can find the Endemic Equilibrium Point? in these models for Covid-19? Or they can be used for data fitting? For my model SVEIRS Covid-19 model if I want to use the A (asymptomatic infection) but our date missing this category, could we assume this category as zero in my data fitting? Is it better to assume the denominator in the force of infection to be (N[t]-D[t]) than the current one (S+E+I+A+R+V)? The wolfram data repository is missing the country "United Arab Emirates"? could you recommend any other data source? Regards

The endemic equilibria can be found with Solve, the way you have been trying. I may just take a long time to get a fully symbolic solution. Unless you want to prove that an endemic equilibrium is unstable when R0 < 1, you generally don't need to determine it symbolically.

For data fitting, it's easier to do things numerically, and an endemic equilibrium can be included in the fitting data.

Regarding the asymptomatic compartment A, if counts of asymptomatic people are not available yo just don't include that compartment in the fitting data. Assuming 0 and would be a mistake.

When using standard incidence for the force of infection, the death compartment(s) are excluded from the denominator, just as you suggest. You can do as N[t]-D[t] or as S[t]+E[t]+A[t]+I[t]+R[t]+V[t].

I'll have to poke around for a data source for UAE cases and deaths.

Posted 1 month ago

Greetings Mr. Jeremy, Greetings... in session2 notebook, how can the "Graph formatting" code at the end of the notebook be modified to include the natural death rate from each compartment, and the death rate from the infected compartment also? (infected compartment includes t outflow: the natural death rate and disease mortality rate.

gstyle[str_] := Style[str, Black, 18]
makeEpiModelGraph[info:{Rule[_DirectedEdge, _String]..}] := Block[
    {edges, g},
    edges = info[[All,1]];
    g = imakeEpiModelGraph[info, edges, Left];
    If[!Less@@Subtract@@(MinMax/@Transpose@(VertexCoordinates /. AbsoluteOptions[g, VertexCoordinates])),
       imakeEpiModelGraph[info, edges, Top],
imakeEpiModelGraph[info_:None, edges:_DirectedEdge|{_DirectedEdge..}, orientation_:Top] := Graph[
       VertexShapeFunction-> Function[{xy, v, wh}, Inset[Framed[gstyle[v], Background->LightBlue], xy]],
       EdgeLabels -> (info /. Rule[edg_, str_String] :> Rule[edg, Placed[gstyle@str, .5]]),

Thanks and regards

Hi Manal, sorry for the delay in getting back to you. While you could go about modifying the code I provided for formatting the graphs in the first few lectures, I think I would recommend making use of the functions in Bob's package and use CompartmentModelGraph instead if you want to be able to automatically include things like death rates to your graph. It would be reinventing the wheel to certain extent to modify my code when Bob has already handled it!

Hi Mr. Jeremy, Thanks for your feedback. Really I did via Dr. Robert's package, but it seems something went wrong or I missed some critical point if you could guide me? This is the screenshot. Maybe because I used the vital demographic as well? Regards.


Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract