solving ordinary differential equations modeling the loss of immunity is very similar to modeling recovery from the diseasethe amount of people losing immunity is proportional to the amount of recovered patients and the length of the time interval we can therefore write the loss in the category as tr in time twhere is the typical time it takes to lose immunity the loss in tis gain in tthe "budgetsfor the categories therefore become nc vts trn ( nc vts ti ( nc ti tr ( dividing by and letting gives the differential equation system vsi ri vsi ir ( ( ( this system can be solved by the same methods as we demonstrated for the original sir model only one modification in the program is necessaryadding nu* [nto the [ + update and subtracting the same quantity in the [ + updatefor in range(n_t) [ + [ndt*beta* [ ]* [ndt*nu* [ni[ + [ndt*beta* [ ]* [ndt*gamma* [nr[ + [ndt*gamma* [ndt*nu* [nthe modified code is found in the file sir py setting to daysreducing by factor of compared to the previous example ( : )and simulating for days gives an oscillatory behavior in the categoriesas depicted in fig it is easy now to play around and study how the parameters affect the spreading of the disease for examplemaking the disease slightly more effective (increase to and increasing the average time to loss of immunity to days lead to other oscillationssee fig incorporating vaccination we can extend the model to also include vaccination to this endit can be useful to track those who are vaccinated and those who are not sowe introduce fourth categoryvfor those who have taken successful vaccination furthermorewe assume that in time interval ta fraction pt of the category is subject to successful vaccination this means that in the time tpts people leave from the to the category since the vaccinated ones cannot get the diseasethere is no impact on the or categories we can visualize the categoriesand the movement between themas