# Model of Epidemic Kinetics with a Source on the Example of Moscow

Authors: Andrey V. Borovsky; Andrey L. Galkin
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

A new theoretical model of epidemic kinetics is considered, which uses elements of the physical model of the kinetics of the atomic level populations of an active laser medium as follows: a description of states and their populations, transition rates between states, an integral operator, and a source of influence. It is shown that to describe a long-term epidemic, it is necessary to use the concept of the source of infection. With a model constant source of infection, the epidemic, in terms of the number of actively infected people, goes to a stationary regime, which does not depend on the population size and the characteristics of quarantine measures. Statistics for Moscow daily increase in infected is used to determine the real source of infection. An interpretation of the waves generated by the source is given. It is shown that more accurate statistics of excess mortality can only be used to clarify the frequency rate of mortality of the epidemic, but not to determine the source of infection.


---

## Body

## 1. Introduction

The situation with the epidemic of the coronavirus infection dictates the use of laser physics methods to describe the development of the epidemic. Elements of a physical model of the kinetics of atomic level populations [1] for active media of laser systems can be used to construct an epidemic model. The four elements are as follows: states and their populations, rates of transitions between states, a time integral operator, and a source of influence. As applied to the description of the epidemic, the population of the settlement (a set of atoms) is divided into groups (atoms are in different states) as follows: healthy nonimmune, infected, recovered, and fatal members of the population. Between these states, similarly to laser physics, transitions take place at their own speeds. The concepts of rates of infection, recovery, and transition to the lethal group are introduced. The rates of emptying of the states of the atom are equal to the reciprocal of the decay time of the states; therefore, following a physical analogy, in epidemic kinetics, the rate of recovery is equal to the reciprocal of the time period of the disease, the same with the fatal outcome. The integral operator in the kinetics of populations realizes the spatial connection of atoms; in the model of epidemic kinetics, the integral operator in time takes into account the infection of members of the population during the incubation period of the disease and mortality outside the incubation period. The time integral operator of epidemic kinetics is reduced to lagging terms, which distinguishes it from the classical SIR model (literally means susceptible-infected-recovered).Epidemic modeling work is divided into two types. The first [3–7] uses a balanced approach applicable for locality or region. In [4], a 7-fraction model of the hepatitis B epidemic is considered. In [5, 6], 6-fraction (SEIARD and SEIAHR) COVID-19 models were studied. In [7], a model with vaccination of the population (SVEAIR) is presented.The second approach [8] further describes the spatial propagation epidemic. Thus, the mathematical model of the Baroyan-Rvachev influenza epidemic [9, 10] is classified as epidemic dynamics, since it can describe the undulating spread of the epidemic between cities when taking into account the passenger traffic (the model is similar to the equations of hydrodynamics). Moreover, all models do not take into account the source of infection [11]. The need to take into account the source of infection is indicated in [2]. Only in this case is it possible to build adequate epidemic models that coincide with real statistics.The aim of this work is to investigate the role of a nonstationary source of infection in delayed epidemic kinetics. In this case, a mathematical model of the epidemic is applied, using elements of physical models of laser physics associated with the description of active laser media.


## 2. Numerical Investigations of the Equations of Epidemic Kinetics with Delay with a Constant Source

The model is constructed on the basis of balanced equations, similar to the physical problem of atomic level kinetics.Let us first consider the model of epidemic kinetics without taking into account the process of vaccination of citizens. For example, in Russia, it was possible not to take into account the vaccination of the population until May 2021. The model considers the development over time of an epidemic of viral infection in a locality whose population is a constant valueN0. The population of the locality is divided into 4 categories: N1, healthy members of the population who do not have immunity; N2, infected; N3, acquired immunity (recovered); and N4, members of the population who fell into the lethal group (died) during the entire epidemic.For a closed population, the law of population conservation is fulfilled(1)∑i=14Ni=N0.To fulfill (1), it is necessary that the natural birth rate is equal to the natural mortality of the population, and the speed of departure of the population is equal to the speed of arrival in the locality. For example, the population of the city of Irkutsk has remained at a constant level for 30 years. In addition, the time of epidemic development T~2-3 years is significantly less than the time of migration processes for settlements such as an average city. All the above values depend on time Nit, and N0 is a constant value.The characteristics of the epidemic development are “global statistics variables,” which are called the number of infected members of the population since the beginning of the epidemicN−(t), the number of recovered members of the population since the beginning of the epidemic N+(t), and the number of the lethal group since the beginning of the epidemic Nc(t). It is easy to see that the values N+(t) and Nc(t) coincide with N3t and N4t in the epidemic kinetic model, which does not take into account vaccination.
(2)N+t≡N3t,Nct≡N4t.The variable of global infection statistics is related to the parameters of epidemic kinetics as follows:(3)N1t+N−t=N0.From here, we get(4)N−t=N2t+N3t+N4t.We see that the number of infected persons can be expressed in terms of global statistic variables(5)N2t=N−t−N+t−Nct.Let us write down an expression for the rate of change in the number of infected persons(6)dN2tdt=dN−tdt−dN+tdt−dNctdt=dN−tdt−dN3tdt−dN4tdt.The first term in the right part describes the rate of infection of members of the population; the second, the rate of recovery; and the third, the rate of transition to the lethal group. The rate of infection will be determined by the number of members of the population who are in the latent incubation period of infectionτ and those who carry the infection in an asymptomatic manner, i.e., undetected carriers of infection
(7)dN−tdt=k−∫t−τtdN−t′dt′Kt−t′dt′+γ1∫0t−τdN−t′dt′Kt−t′dt′.Kt−t′− is the core of the integral operator, which describes the decrease in infected individuals over time. The type of model function Kt can be selected as
(8)Kt=exp−tτ+2.The real form of the kernel of the integral operatorKt should be determined by comparison with the statistics of the epidemic.AssumingKt=1, which somewhat overestimates the number of persons transmitting the infection, integrals in (7) are revealed and lagging terms appear
(9)dN−tdt=k−N−t−1−γ1N−t−τ.The rate of recovery of infected citizens is proportional to the number of infected at a given time(10)dN3tdt=k+N2t.The rate of death is determined by the number of infected persons, taking into account the delay in the development of the disease(11)dN4tdt=kcN2t−τ.Given the large volume of testing of the population and the isolation of asymptomatically ill persons in the city of Moscow, it is possible to putγ1≅0. A significant change in the values N3t and N4t occurs at times t>3τ, while N2t is sensitive to changes at times t<τ. Therefore,
(12)N−t−1−γ1N−t−τ≅N−t−N−t−τ≅N2t−N2t−τ.As a result, we come to the following modification of the epidemic kinetics model(13)dN1tdt=−k−N2t−N2t−τ−А,dN2tdt=k−N2t−N2t−τ−k+N2t−kcN2t−τ+А,dN3tdt=k+N2t,dN4tdt=kcN2t−τ.The equations with lagging terms in (13) are not first-order differential equations. The delayed term N2t−τ can be decomposed into a Taylor series. At the same time, higher-order derivatives appear. For example, in [2], the second equation in system (13) was reduced to an ordinary differential equation of the second order, which had a bell-shaped solution of the epidemic wave type.Here, the coefficientsk−,k+,andkc have the dimension of frequency, i.e., the inverse time, which is measured in days. All data in epidemic statistics are presented per day. The number of different categories of the population is determined in persons.The valueA represents the rate of quasistationary permanent sources of infection. The value of A can change slowly. However, the time of such change is significantly longer than the inverse frequencies present in the model and the time of the latent incubation period of infection development
(14)T≫k−−1,k+−1,kc−1,τ.We note the following circumstance, which negates the choice of a nonstationary model for epidemic kinetics. Numerical calculations have shown that the unsteady mode of epidemic development occurs only at the initial stage within 90-100 days. All the subsequent time, the epidemic develops in a quasistationary mode, which weakly depends on the expression for the rate of infection. (Such modes are well known in laser physics. They arise for long or generally stationary pulses of pumping the active medium of the laserAt≈Const.)For the COVID-19 disease caused by the SARS-CoV-2 coronavirus (as opposed to, for example, the flu), we take the recovery time of theτ+=15 Day, then the recovery rate is the reciprocal κ+=1/15Day−1. The frequency of death is similar to κc=1/50Day−1. The duration of the incubation period of theτ=7 Day and the frequency of infection are represented by the expression
(15)κ−t=ptΔṅt=ptΔṅ0tN1tN0,∑i=14Nit=N0.Expression (15) takes into account that the rate of contact of an infected person with uninfected persons should decrease as the epidemic proceeds in proportion to the share of healthy noninfected persons in the population of the locality
(16)Δṅt=Δṅ0tN1tN0.Here,Δṅ0t− is the rate of contacts with all persons who were not ill and recovered.The following expressions are used for the frequency of contacts and the likelihood of infection:(17)Δṅt=Δṅin;0<t<t′Δṅin−Δṅin−Δṅft−t′t′′−t′;t′≤t≤t′′,Δṅf;t′′≤tpt=pin;0<t<t′,pin−pin−pft−t′t′′−t′pf;t′′≤t.;t′≤t≤t′′,Formula (17) describes the smooth introduction of sanitary standards. Designations Δṅin,Δṅf and pin,pfcorrespond to the initial and final values of the rate of contacts and the possibility of infection. The calculations used the valuesΔṅin=30,Δṅf=20 and pin=0.01,pf=0.005.The statement of problem (13) should be supplemented by the condition of nonnegativity of the sought functions and by the initial conditions, for example,
(18)Nit≥0,i=1,..4,N10=N0,N20=1,N30=0,N40=0.To numerically solve the problem of epidemic kinetics (8)–(12), a difference scheme is used, a feature of which is memorizing the prehistory to a given moment t (for example, with a step of 0.1 days) of already obtained functionsNis,s≤t.In the approximationN1t≈N0, the κ−t in (15) does not depend on N1t and the second equation in (13) becomes independent, and the remaining components are determined through this solution N2t.In the approximationN1t≈N0 at small times without taking into account the delay from the second equation of system (13),
(19)dN2dtt=κ−0−κ+N2t+А.We get the analytical expression
(20)N2t=N20ek−0−k+t+Ak−0−k+ek−0−k+t−1.An outbreak of morbidity (20) with a positive incrementk−0−k+>0can be extinguished only due to the lagging terms in (13).In a stationary case,(21)dN2dt=κ−N2t−N2t−τ−κ+N2t−κсN2t−τ+А=0,N2t=N2t−τ.The solution, as in stationary models of laser physics with constant pumping of the active medium, takes the form(22)N2=Ak++kс.The stationary value (22) does not depend on N0,k−.Figure1 shows epidemic curves whenN0=6∙105 a constant source of infection operates in a settlement with a population of people (a medium-sized city). Such a source can be an infected contingent arriving in a settlement with the external passenger traffic. Curve 1 corresponds to the source A=10 people/day and t′=20 and t″=50 in coefficients (17). A later adoption of quarantine measures t′=50 and t″=80 with the same source corresponds to curve 2, and the outbreak is extinguished only due to the lagging terms in (13) before the quarantine measures are taken. Both curves tend to one stationary value (22) determined by the action of the source. Curve 3 corresponds to the source A=50 people/day and t′=20 and t″=50 in the coefficients (17) and reaches a stationary value (22) 5 times greater than curves 1.2. Note that the wavy behavior to the right of the maximum is due to the complex structure of the second equation of system (13). The calculations do not depend on the valueN1; this dependence can manifest itself when the entire population of the city is ill or when it is vaccinated.Figure 1
Solution of system (1) withN0=6∙105: A=10, t′=20, and t″=50 (1); A=10, t′=50, and t″=80 (2); A=50, t′=20, and t″=50 (3).The form of the morbidity curve is standard (rise-decline), qualitatively coincides with the results [5, 10, 11], and differs significantly over long periods: an asymptotic tail instead of a sharp exit to 0.


## 3. Inverse Problem of Epidemic Kinetics

To analyze the properties of the solution to problem (8)–(12), it is necessary to know the function of the real source At. This function can be determined by solving the inverse problem of epidemic kinetics. In laser physics, this method is used to determine the pump function from the experimental emission spectra of the atoms of the medium. In epidemic statistics, infection rates (per dayΔt=1) and death rates (per day) are used, which correspond to the solution components (13).
(23)ΔN2Δt=κ−N2t−N2t−τ+А,(24)ΔN4Δt=κсN2t−τ.Determination of the source of infectionAt can be done by solving the inverse problem (8), where the left side is the statistics data [12] and the right side is the solution to system (1) with the desired sourceAt. In the formulation of a direct problem, the equations and coefficients are unchanged, and only the initial population size changes. Step-by-step time fitting is used as a numerical method for solving the inverse problem. Figure 2 shows the daily statistics of morbidity in Moscow (N0=1.2∙107) for the 539 days starting from 12 March 2020 (1, solid line); the source determined by solving the inverse problem (23) and (13) (2, dashed line); calculated data (23) of model (1) with the obtained source (3, points).Figure 2
Statistical data of daily morbidity in Moscow for 539 days starting from 12 March 2020 (1, solid line); the source determined by solving the inverse problem (23) and (13) (2, dashed line); calculated data of the model (1) with the obtained source (3, points).The average deviation of the calculation data and statistics is 629 people/day, which corresponds to a relative error of less than 7%. A wave can be understood as a maximum amplitude with a width determined at half height. Waves in the number of infected, the daily increment of the infected, and the source of infection have shifted maximum positions. The structure with source waves is primary and determined according to epidemic statistics. Two main factors for the emergence of waves of infection are as follows: the importation of infected with an external transport stream and an internal increase in contacts due to noncompliance with sanitary standards. The combination of these factors is indirectly taken into account in the source as a result of solving the inverse problem. The first wave has a maximum of 54 days of calculations (5 May 2020), a rather narrow width due to the adoption of strict quarantine measures in the form of lockdown, which are taken into account in the quasistationary mode in the source indirectly through statistical data. The second factor seems to have been predominant. The second wave, which was more extended, had a large amplitude and was formed by a combination of two factors. The development of the third wave is associated with the emergence of the Indian strain, and the decline of the wave is due to vaccination.Figure3 shows the curves of morbidityN2t (solid) and lethalityN4t (dotted line), calculated according to model (13) with the obtained source and the modified coefficientkс=1/200.Figure 3
Estimated curves of morbidity (solid) and mortality (dotted line) with the received source and the modified coefficient.The integral number of infected people in Moscow since the beginning of the epidemic of 1 September 2021 is the following: statistical data, 1,568,767 7 people; calculation data from the received source, 1,597,000 people.Statistical data (23) and (24) in the quasistationary approximation taking into account (22) can be represented as
(25)ΔN2Δt=κ−N2t−N2t−τ+А≈А,(26)ΔN4Δt=κсN2t−τ≈kсAk++kс.The quasistationary approximation is well known in the theory of laser media. Such a regime arises when the time of action of the pumping of the laser medium is many times longer than the times of depletion of the energy states of the atoms of the active medium.Both statistics (25) and (26) are determined by the source of infection A. The model gives overestimated values for the lethal group, which requires correction of the coefficientκс. The coefficient κсis reduced by 4 times, while the residual value changes from 626 people/day to 644 people/day with the same source of infection, i.e., the dependence of the model results on the coefficient is not significant. Mortality statistics (26) cannot be used to independently determine the source of infection due to the uncertainty in the frequency of deathκс.


## 4. Analysis of the Correlation between the Source of Infection and the Number of InfectedN2t

The sourceАt is obtained as a result of solving the inverse problem of epidkinetics and the N2t solution of the direct problem of epidkinetics with the obtained source. Figure 4 shows a graph of the parametric АN2as follows: thick curve, the first wave; thin curve, the second wave; the dotted line, the third wave; curve points–stationary relation
(27)A=N2k++kс,kс=1200.Figure 4
Parametric graphАN2: thick curve, first wave; thin curve, the second wave; the dotted line, the third wave; points–stationary relation (9) with a coefficientkс=1/200.The dependence is of a spiral nature, so that the valuesА,N2can be interpreted as waves, but their maxima do not coincide, i.e., the waves are different. In general, the values А,N2 are correlated, tend to a stationary connection, with the exception of the lower part of the graphs, where N2decreases at constant A, which does not contradict the course of the incidence curve (Figure 1). CorrelationА,N2means that the source is formed by the number of infected. Areas marked with numbers 1, 2, 3a, and 3b are highlighted, in which the source sharply increases with the number of infected. Section 1 is a jump on the 50th day of observations (1 May 2020, a holiday in the Russian Federation). Plot 2 - smoother growth 171-202 days of observations (30 August 2020 – 30 September 2020 - the end of summer vacations and the beginning of the school year). Points 3a -393 (12 April 2021), 3b - 456 (11 June 2021). A jump to the holiday, the second factor works, 2 is a combination of factors, 3a is the second factor, and 3b is the first factor (the import of an Indian strain).Application of the theory for other cities gives [13] for 1 September 2021 the position of the coordinate of the end point of the spiral: Berlin, A=906; New York, A=4126; and Catalonia, A=636. Quantity information N2is not available. Recession of the epidemic is the approach of the spiral to the origin.


## 5. Accounting for Vaccination in the Model of Epidemic Kinetics

Vaccination above is indirectly accounted for by reducing the source of infection. Direct accounting of vaccination can be performed within the framework of the following model:(28)∑i=14Nit=N0,dN1dtt=−κ−N2t−N2t−τ−Аt−Bt,dN2dtt=κ−N2t−N2t−τ−κ+N2t−κсN2t−τ+Аt,dN3dtt=κ+N2t+Bt,dN4dtt=κсN2t−τ.This adds the vaccination rateBt with a minus sign to the equation forN1 and a plus sign to the equation forN3. In the stationary case, in the approximationАt≈const,Bt≈const(29)N1=N0−A+Bt;N2=Ak++kс;N3=k+Ak++kс+Bt;N4=kсAk++kсt.Vaccination with this approach reduces the size of group 1, increases the number of group 3, and directly not affects the abundance and groups 2 and 4. The indirect impact of vaccination on the number of groups is obtained from (29) in the form of evaluations and infection rate with the vaccinations.
(30)κ−~N1tN0~1−A+BtN0.


## 6. Reproduction Coefficient

In epidemiology, the process of epidemic development is often described using the reproduction coefficient. The rate of spread (reproduction) of virusRt is an indicator that determines the average number of people infected by one patient before his isolation (or recovery). It is calculated based on the data on the increase in new cases over the last τ0 days. The virus spread coefficient is used to make decisions about the transition to the next stage of lifting restrictions
(31)Rt=1N2tdN−tdt∙τ0.N2t is the number of infected persons, dN−t/dt is the rate of infections per day t, τ0 is the time of the spread of infection by an infected person, and N−t is the global statistic parameter (the number of infections since the outbreak). Figure 5 shows the time course of the reproduction coefficient based on statistical data for the city of Moscow and the value of τ0=τ+=15 days.Figure 5
The time course of the reproduction coefficient (31) for the city of Moscow. Time t is specified in days.The graph indicates an obvious weakening of the viral infection compared to the first wave of the virus, as well as the transition of the epidemic to a quasistationary mode, where(32)Rt≈τ+k++kс≈1.


## 7. Discussion

The study of epidemic kinetics within the framework of a new theoretical model reproduces the course of the infection curve with three characteristic areas: rise, fall, and exit to a stationary value. The spread of infection is described by paired interactions with average speeds and the source of infection. The dynamic system underlying the model is nonstandard—with lagging terms and a source that is determined as a result of solving the inverse problem. The presence of a source leading to the possibility of describing the stationary mode of the epidemic distinguishes the proposed model from all other previously known ones.The model with a source is minimal in the number of group states, which allows, with the uncertainty of the coefficients of the problem, to set and solve the inverse problem of epidemic kinetics to determine the time course of the source of infection for Moscow by one set of statistics in the form of a daily increase in morbidity. Another set of statistics is used to refine the coefficients. Next, a direct problem is solved with the found source and all the characteristics of the temporal evolution of the epidemic are obtained, including an estimate of the reproduction coefficient. Solving the direct problem of epidemic kinetics with the found source, the authors obtained a coincidence with the statistical data for the city of Moscow with an error of not more than 7%. The formulation and solution of the inverse problem of epidemic kinetics within the framework of the model under consideration is the second distinctive feature of the approach used. The first wave of the epidemic proceeds in a nonstationary mode and has the duration of 90-100 days. Further, while maintaining sanitary measures, the epidemic proceeds in a quasistationary mode, for which the source of infection is crucial. The spiral relationship of the source with the number of infected determines the dependence of the source onN2t and the wave-like nature of the epidemic development. The beginning of the rise of waves is associated with an increase in contacts and incoming traffic flow, and the decline of waves is associated with the introduction of quarantine measures. The impact of vaccination in Moscow on the reduction of the third wave is estimated as follows: mild quarantine measures reduce the level of infection by 3 times, vaccination gives an additional lowering coefficient (30), where the level of vaccination reaches values of B≈105person/day.If appropriate statistical data are available, the model can be implemented for any region of the world. The data given for Berlin,A=906; New York, A=4126; and Catalonia, A=636 are not complete in this sense. They show only one coordinate on the spiral in space N2,Аat a fixed moment. The spiral connection through passenger traffic extends to other cities and regions. Individual turns of the spiral may be degenerate due to the small amplitude of the wave. The ultimate goal is to get out of the spiral to the origin: N2=0,A=0.Note that our model (13)–(18) is descriptive. In order to proceed to predictions, you need to study in detail the formation of the structure of the source At. We hope that this part of the problem will be solved in future works.


## 8. Conclusions

The article presents a study of epidemic kinetics within the framework of a new theoretical model using kinetic elements for populations of atomic levels of an active laser medium as follows: a description of states and rates of transitions between them, an integral operator, and a source of influence on the system. In this model, the time integral operator is reduced to an algebraic form with lagging terms. It is shown that the concept of the source of infection should be used to describe a long-term epidemic. With a constant source of infection, the epidemic, in terms of the number of actively infected people, goes into a stationary mode, which does not depend on the population and the specifics of quarantine measures.Structure of the source of quasistationary infection is not known; however, the task of describing infection statistics for a locality (city) can be solved as follows. At the first stage, the inverse problem of epidemic kinetics is solved. To determine the real source of infection, statistics were used in Moscow on the daily increase in the number of infected. At the second stage, the direct problem of epidemic kinetics with the already found source is solved. This made it possible to obtain epidemic curves that coincide with statistical data with an accuracy of 7% for the entire long-term stage of the epidemic development, taking into account three waves of infection.The interpretation of epidemic waves generated by the source is given. Waves of infection are formed by two main factors: external incoming passenger traffic and, possibly, an increase in contacts between citizens due to the weakening of sanitary measures. In the formation of the first and third waves, an increase in contacts prevails, and in the formation of the second, a combination of factors.It is shown that in the quasistationary mode of epidemic development, the use of two statistical data on the frequency of infection and mortality makes it possible to determine the epidemic mortality rate, which for Moscow was about 0.5% per day of the number of infected. At the same time, accurate mortality statistics cannot be used to determine the source of infection due to uncertainty in theκс mortality rate.



---

*Source: 6145242-2022-02-18.xml*