Thanks for adding the additional detail about the model, it helps one understand the results better.

It looks like you are using the data from Hubei, based on the large number of confirmed cases. However, when I downloaded the JHU data from github (https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data/csse_covid_19_time_series) last evening, I got similar but noticeably different values. What is your data source, and what preprocessing did you do? What was the adjustment that you made?

What assumptions did you make when choosing a susceptible population size of 180,000? The population of Hubei is very much larger:

In[22]:= WolframAlpha["population Hubei China", "Result"]

Out[22]= Quantity[58160000, "People"]

It will have a dramatic effect on the dynamics, especially if one uses mass action incidence for the force of infection as opposed to standard incidence. Which did you use? Models for the SARS epidemic used mass action incidence.

I have found these data quite challenging to model. My current model has 7 compartments and includes quarantine. It's still very preliminary, and I hope to have it ready next week, but here is a peek at its structure (as a Petri net);

I made a notebook deriving the logistic model from your SEIR equations with a simple assumption.

Thank you, your points are well taken and I will do so in the next round ... there are already some updates and some changes and some things that I still need to clarify (or else I am creating confusion) ... but I think I will wait a few days and see how close future data is to the model without having modified it further ... I am hoping in the end to have an accurate description of what is observed by the data and I hope to publish the equations that give the model. Thank you once again.

You are of course right about all of what you say. My goal was to try to find a phenomenological reason with implications about the values of the parameters that would make it possible for me to model the data, and I succeeded in that sense. I agree with you in that one has to be cautious about the scaling. But for the data that we have, we have some kind of model. I will update this data daily and see if it continues to fit the model. I expect it will and what the model predicts is in line with the expectations of some epidemiologists.

Right, real data is messy ... you have to do the best you can with it. I think that what I did somewhat works.

About recovery, it should mean "not infectious" given the tests that the patients need to go through to be discharged. However, there are cases documented when discharged individuals test positive again later on, leading to the suspicion that the disease might be biphasic in some instances.

And correct, each focus will probably need its own, if similar, parameters. I am tracking several of the new foci, and we will see what happens with time. Unfortunately, it is harder to get data broken up into clusters of interest. Surely it can be obtained, but it is not easy to get it.

Hello,

Let me try to explain how I tried to model these data and try to address some of your questions, not necessarily in the order they appear. First, I decided to work with the simplest model possible, to keep the number of parameters as low as possible. Otherwise, estimating the parameters gets quite difficult. It seemed to me that the basic SEIR model should work, if only we modified it to reflect the situation at hand. Then I made some thought experiments to try to make these modifications and parameter estimations. Let me explain a couple of these with an example.

The first thing I tried to understand is the effect of the containment measures enacted on the parameters of the basic SEIR model. I illustrate what happens with an simple example. Suppose you have ten families with four members each and that you all of a sudden detect three infected individuals, each in a different family. Now suppose you then immediately put each family in a different house, and you prevent them from having any contact with each other. The effect of this is essentially to remove 28 individuals from the susceptible group, which would otherwise be 40. Only 12 individuals, those families with an infected member, are now susceptible. Understanding this told me that I might have a chance of modeling the data if I lowered the susceptible population to an adequate number. I then used further information about the other possible value of the parameters, including information about maximum observed incubation periods. And then it was a matter of some fine tuning to get the right fit for the I curve. I observed that even if the actual data was only a sample of the real number of infections, the shape of the curve should be the same with the real numbers. You just multiply by a constant factor. So say that, as suggested, there might have been ten times (or twenty) times the number of cases as reported, then the susceptible population would have been 1800000 or 3600000 and likewise, you would have ten times more infections, recoveries, etc. The key is getting the shape of the curve right.

In the model, I tried to get the I curve right, and then, modify the equation for R with a root factor to account for a delay of the onset of recovery that fits the recovery data (it takes a while for people to recover).

Now about the data. I started early on, and I only had data that I could copy by hand ... and then I just kept doing it that way, among other reasons, because it forces me to see what is going on, rather than just letting the computer read it. From the 13th of February on, I subtract a constant of about 14000 to the JHU data, because of the change of counting method that lasted only a few days and the leap in the numbers. This is the simplest quick and dirty fix that keeps the data in line with the trend and nullifies the leap in the count which occurred on that particular day.

I am quite confident that the model is basically working well. As the days go by, I will continue to check it and adjust it.

If you think about the most difficult issue regarding action to contain the spread, given that you can't test extensively, is to decide when you have a large enough number of detected cases to warrant some containment action. It is a difficult decision to make, and the better data and more extensive testing you have, the better prepared you are to make a wise decision. It would be very useful to have a very good estimate of what the real number of cases is given a number of detected cases that you might have as a sample. We don't really know. Some literature has suggested that it might be 20 times the number of detected cases.

It might be that in the new foci of infection, the model might have to be modified to fit the circumstances of the situation. I am following several foci. I did model the cruise ship situation, and the with almost the same parameters I managed to obtain numbers that were quite very close to the actual numbers while that outbreak lasted. When I have time I will revisit that model and see if the improvements in the model here can give even better results. I need to consider all days, I don't have all the days.

I hope this answers some of your questions. I hope I get around to cleaning the notebook and providing all the details ... but in a nutshell, it is the simplest possible SEIR model you can think of, almost.

Hello,

I will post this soon for the SIR models.

Hi Enrique,

Could you possibly let me know how to apply and lift restrictions? Many thanks, Martin

Look at the equations AND the initial conditions (for S). What you see (relative to the y axis) is what you get ... (I should include a plot with S) ... hope this helps; if not, we can come back to it at some point.

Thank you!

I would like to add that in China the situation was strictly controlled until the end, so what had to stay constant (susceptibility) did ... I am afraid that elsewhere this will not happen. I already have a model which allows for the growth of susceptibility. It is beyond the scope to discuss it here, but I plan to make a notebook available in the aftermath, when we have the whole picture. So when restrictions are lifted, susceptibility, the parameter that controls the size of the epidemic, grows again, maybe ever so slowly, all depending ... all the measures put together work to determine that one quantity, which is the overall most important quantity.

Enrique

Enrique,

Attached are two notebooks which I am currently using. The fitting methods are embedded in the report procedures. When fitting to cumulative cases, the large numbers in the latest data dominate the fit, but they are probably also more accurate. At the end of the CDCData.nb there is a method for fitting the first differences of the cumulative data to a first difference formula derived from the logistic model, if you also want to see more directly the influence of the early data. We are about to enter the tail phase of the epidemic in the U.S. The model requires an exponential decay in the tail. That didn't happen with South Korea, but the China data did fit. If the U.S. data doesn't follow the tail behavior, that would indicate the model isn't valid.

Bob

Attachments:

CDCData.nb

KeyStates.nb

In this section I discuss an SIR-like model which works almost as well as the SEIR model, at least with the Chinese data. I will have to investigate whether I can make it work as well with other data sets.

Our SIR like model is described as follows

s'(t) = - Beta * s(t) * i(t) / p,

i'(t) = Beta * s(t) * i(t) / p - Gamma * i(t - n),

r'(t) = Gamma * i(t - n)

The function s(t) is the number of susceptible people (the people that can get exposed to the pathogen) at time t; i(t) is the number of people who are infected and infective; r(t) is the number of people who have become resistant to the pathogen: they have recovered and developed immunity or died. Now the parameters.

beta is the rate of infection or "force of infection", gamma is the removal rate, and n is a shift parameter used in part to line up the curves (see a bit of an explanation of this in the SEIR section). The parameter values are in the titles of the pictures for each country. In general we assume i(0)=1 unless stated otherwise in the model label. Also, s(0)=p, and r(0)=0, unless otherwise stated.

We present in the picture, an SIR model and its parameters that fits the Chinese infection data. With time, I will try to fit the SIR-like model to other data sets. I also plan to investigate the known analytic solutions (although I need to understand how the shift parameter changes the classical solutions) in order to attempt automatic fitting via computational optimization in Mathematica. This is work in progress.

August 30: Positivity rates picture updated. For the time being, we will no longer update this picture

August 10, 17: Positivity rates picture updated. Mexico, which has a huge positivity rate, aside, Sweden still has a very high one.

June 15,22,29: Updated June 29, July 13, 27. The positivity ratios picture will not be updated again until August 17, or then again, only occasionally. All countries in the picture have brought this ratio down over time, except Mexico.

June 1,8: Positivity rates updated

May 21: The positivity rates file will be updated on Mondays from now on. It is updated today.

April 21-May 12: The SIR models will now only be updated occasionally. The positivity rates picture will be updated daily.

April 19-20: Finland updated. A new SIR model for Finland with a different recovery schedule is included. We include now the current positivity rates for various countries (number of cases/number of tests)

April 18. There is now a SIR model for Finland which uses JHU data

April 14: I have added the SIR model for Italy. I will not update this daily. Later, I will make a notebook available for this model. If I am able to complete it, I will also have in the notebook a program to find the parameters to fit the model to the data - but I cannot promise I will have that, at least not soon. Note the forecast is somewhat more optimistic, both in regard to total susceptibility and duration of o

Attachments:

FiSIRmodel6bfr.png

FiSIRmodel6c.png

ISIRmodel6.png

SIRmodel6.png

SPositivityRatesToday.png

In this section I attach the notebook. I will update it soon and add new features to it, such as the effect of lifting restrictions early.

August 30, September 6: the pdf document has been updated. Next update at the end of September

June 29, July 5, 13,19, 27, August 2, 10, 17: the pdf is updated. Next update in a week

This notebook is also posted in a reply to Kaurov, above the Scandinavian countries section. I will post (hopefully soon, I did not get it ready June 22) the new notebook for adjusting parameters automatically in both sections. In the reply to Kaurov, there is also a pdf document with daily new cases curves for several countries, now also attached here. This will not be updated until August 17, or then again, maybe only occasionally. It is interesting to see that some countries which brought their numbers down initially have settled in a plateau of a small number of daily cases, without being able to completely eliminate the v

Attachments:

SEIR model C 27 ...nb

SmoothDaily.pdf

I'm afraid I also lacked that foresight and only have a little saved (and most of that's just for the USA) as it was a while before I realized the significance that the changing variables had. I've begun a more thorough archive NOW though and I expect to keep it up for the future! I'll be happy to hand that archive to you whenever you like.

It strikes me that the hand-adjustment of the variables that you've been doing is analogous to the best-fit line trick of laying a ruler through a scatterplot of data with a linear trend and eyeballing and equalizing dots above/below the line. It is quite effective for datasets with a clear linear trend and likewise, your methods are self-evidently quite effective for these SEIR-like trend datasets.

A much more computationally intensive method would be to apply the appropriate Euler's method (or Runge-Kutta, I think 4th order? method) to the SEIR equations for the dataset at some early time per each group reporting data (world/country/province) with enough volume and shape to be fit to, and re-do that best fit for each day data is reported. The changes required of the variables to maintain that best-fit would then be tied to a more mathematically rigorous definition of "best" fit as well. Doing this and comparing the changes of variables over time per group would also seem likely to help determine correlations between generalizations of how the epidemic progresses as well as comparisons of each group's susceptibilities and methods of combating the virus.

Doctor Moreno, I've been following this discussion, which I found through Robert Nassar's, for nearly a month and I've been consistently impressed by not only your research and models (which are the best fits to the data that I've seen ANYWHERE so far) but also by the explanation and discussion of the work. I actually made a Wolfram account primarily so that I could post this message to let you know how much this kind of work (and education about the underpinnings of the work) is appreciated. I look forward to further exploration of Wolfram's community in hopes of finding similar quality discussions of subjects more closely related to my field, as well as to seeing the full body of this work in publication. Many thanks!
As I understand it, you're plotting a standard SEIR curve and doing a best fit to the reported data (as reported by various health organizations)- therefore looking ahead and using this model predictively the future curve is an indication not necessarily of how many are actually infected but of how many will be reported as such on each coming day. To me this is very interesting because it sidesteps many of the underlying issues of tainted data due to under testing, selective testing, flaws in tests, etc. and just rolls them up into the variables of the model. Picking those apart after the fact is very interesting, but in the moment, given the data set(s) available, this seems like perhaps the most useful utilitarian presentation of a model.
It would be very interesting to me to see as a 3D plot, or even as an animation (use time to show another variable?) how your model has been adjusted as time and the epidemic has progressed. I've done that crudely for myself just by flipping through the images I've saved for reference of the plots you've presented and updated daily watching the curves change as new data points appeared and I think a representation of that with the variable values adjusting as the model adjusts would be revealing.
I think of most things as plots of curves or surfaces first, equations driving their shapes second so I tend to be biased toward a graphical representation. Therefore, I'm highly interested in a plot of the underlying variables used for each best fit SEIR curve day by day as the model was adjusted to fit each new data point. Having a near-constant value the whole way through would be the classic case, where the model is an excellent fit to the phenomenon all the way through. In this case though, I believe that one of the big differences is in your use of a much smaller S initially than the total population of each region reporting data (as one example variable). The explanation that the entire population isn't actually susceptible early on (for instance, due to social distancing/quarantining attempts) rings true to me (I think back to chemical rate or mixing equations from Calc II and remember that those always assumed "well mixed" tanks of solutions, which is NOT always the case in real life!). Therefore, a plot of how the S needed to change with time within each country to keep a best-fit would likely be very revealing of how the epidemic surged (or didn't) through the population and changes in it's slope might correlate with changes in testing, reporting of data, treatment methods, institution of distancing measures, etc.
A deep dive into the other variables' change over time likely would be equally revealing. I'm too new to Wolfram software to know if this "historical" data for earlier fits is still included in this notebook, but it would of interest!

Hello again,

Thanks, give me a day to work on this … I am about ready to quit for the day, except for minor things.

To answer your questions:

1) If you go to the top of the post, I posted the equations and explain the parameters. Beta is the rate of infections, as it is known. Sigma is the rate at which the exposed become infective. Gamma is the removal rate. p is the total number of susceptible individuals (this is equal to the total number of cases there will be at the end). The shifts are used to align the curves with precision with the data and to account for the fact that in different countries, the "release policy" of sick patients varies … I won't go into details, maybe some other time; mathematically, it works.

2) The main effect of the lockdown is reflected in the size of p, as I explained in the discussion with Robert Nachbar above. If you read that whole discussion, at some point I explain why this is the main effect of the lockdown; it is the secret to making this work, and it also makes sense from the epidemiological point of view, not just mathematically (since you are a physicist you might appreciate the thought experiment discussed above)

3) Initially I had a much more complicated model which I used to estimate the parameters for the basic SEIR model more or less automatically, along with the assumption that p had to be low with a lockdown. Then I discarded that and learned to adjust the parameters as needed for the initial data set (China) and then for other data sets. As I said, I made a thought experiment that pointed the way and taught me that the key parameter to make this work is p. If you read the discussion with Robert Nachbar I explain this there. The main effect of the lockdown is to lower p.

I usually fit first when there is already enough data, and modify for a while until the model is stable; then it just works well from a time on and I don't have to adjust the parameters so much, or even at all. I have worked really hard to learn to fit the parameters. If I had time, I would use machine learning to do it, but right now, I have had to learn to do it myself. You could also do it computationally, but the results are not always the best, and that also takes time and work to get it right. I have learned to modify the model for each data set fairly effectively, each data set has its own strangeness, but the data sets are not so different that some basic values for the parameters could not account for it for the most part if you know how to adjust them properly … that is where ML would be useful, and what I have learned to do without ML.

4) In this model you cannot split recovered from fatal.

5) Yes, I will be glad to … give me a day and I will write to you. I am working 24/7 on this and many other related things. So I will get in touch with you tomorrow.

Best, Enrique

Enrique, I am writing from Italy. I am a physicis. Today I was trying to see how to adapt SEIR to the OBSERVED value and I found your work ! Impressive. Congrats. You know how bat the situation is in Italy. Your work would be extremely useful. I would like to suggest to additionally test your work on a couple or regions in Italy: Lombardia, Emilia Romagna, Veneto and Piemonte. In this way you could spot if there are or not differences in the fitting parameters. Would it be possible at all ? Please reply me asap: time is of essence. With my best regards Roberto Battiston

Hi,

I have updated my notebook (March 20) with the China model and data. I will add the other models when they stabilize.

Hello,

The published notebook is not my working notebook. I am on travel for a week. I hope to get to updating the published notebook then, apologies, and thank you.

Thanks!

Dear @Enrique Garcia Moreno E. thank you for sharing! Could you please share the notebook with Wolfram Language code? You can attach it or embedded into the post.

Thank you, it's much more clear now.

As a teacher of graphing and visualization, might I make a proposal that you properly label chart axes with units and also label charted values (which colour is which variable).

This makes it much easier to instantly understand the graph.

Thanks for your efforts and sharing them.