key: cord-239632-dyvhzvmm
authors: McGuigan, Michael
title: Pandemic modeling and the renormalization group equations: Effect of contact matrices, fixed points and nonspecific vaccine waning
date: 2020-08-03
journal: nan
DOI: nan
sha: 
doc_id: 239632
cord_uid: dyvhzvmm

In this paper we find common features between the equations that are used for pandemic or epidemic modeling and the renormalization group equations that are used in high energy physics. Some of these features include the relation of contact matrices in pandemic modeling and operator mixing in the renormalization group equations. Another common feature are the use of flow diagrams and the study of fixed points both in pandemic modeling and in evolution under renormalization group equations. We illustrate these relations through the study of some cases of interest to the current COVID-19 pandemic. These include pandemic modeling with mixing between different age groups and also contact matrices associated with contact between countries. For the final example we study the effect on mortality of waning from nonspecific vaccines which are designed to combat different pathogens but nevertheless may lessen the severity and mortality of COVID-19 infections.

It is well known that in physics phenomena that are governed by the same type of equations can benefit from the sharing of techniques for the solution and visualization of the solutions to those equations. For example the simple harmonic oscillator and electronic circuits are described by similar equations and insight from one equation can be applied to the other. The first order nature of the differential equations used in compartment models for epidemic modeling draw a striking resemblance to the renormalization group equations used to describe high energy physics models. Essentially each compartment or the number of class of individuals in the model correspond to a coupling constant in the high energy physics model. One can even go beyond this and draw a closer relation to scaling structure of virus spreading and evolution of a physics model under change in scale [1] [2] . In this paper we illustrate the relation between. the techniques used in high energy physics and apply these to cases of interest to pandemics with similar characteristics to the current COVID-19 pandemic.

This paper is organized as follows. In section one we give a basic introduction to our approach. In section two we give the basic equations that are used in pandemic and epidemic modeling based on classifying individuals in different compartments. In section three we give the basic renormalization group equations of high energy physics associated with coupling constants. We show how basic solution techniques and visualizations involving flow diagrams can also be applied to compartment models of pandemic modeling and how these can be used to identify fixed points of final state equilibrium in pandemic modeling. We also draw a parallel between the between pandemic contact matrices and operator mixing where a coefficient of an operator that may be zero at some scale nevertheless can become nonzero at a different scale through mixing with other operators in the theory. In section four we illustrate the above relations using contact matrices between different age groups using demographics and pandemic characteristics similar to the current COVID-19 pandemic. In section five we further illustrate the techniques using contact matrices between different countries which encounter the virus through delayed interaction and visualizing the curves for infection and mortality during the pandemic. In section six we study the effect of nonspecific vaccines that do not prevent infection but may lessen the severity and mortality of the pandemic by studying the model with separate age groups both with and without the use of a nonspecific vaccine. Finally in section seven we state the main conclusions of the paper and opportunities for further exploration of the common features of pandemic modeling and remormalization group equations. 

Here S is the fraction of the population of susceptible individuals, I is the fraction of infected individuals and R is the fraction of recovered individuals. In these equations β and γ are constants associated with the epidemics that determine the likelihood of infection and the rate of recovery. Another model that we will use is the SIRD model which is similar to the SIR model but takes into account mortality associated with the infection. For the SIRD model the compartment equations are:

with initial conditions

Here S is the fraction of the population of susceptible individuals, I is the fraction of infected individuals, R is the fraction of recovered individuals and D is the fraction of deceased individuals. In these equations β, γ and µ are constants associated with the epidemics that determine the likelihood of infection, the rate of recovery and the rate of dying from the infection. For small number of compartments these epidemic equations are not difficult to solve numerically using software such as Mathematica or MATLAB [8] . In figure 1 we show the numerical solution for an epidemic described by an SIR model for parameters (β = 3/2, γ = 2/3) with initial conditions (S 0 = .999, I 0 = .001, R 0 = 0). 

In the compartmental models of epidemics one has different compartments associated with susceptible, infected , recovered and deceased individuals. In the renormalization group equations of high energy physics instead of compartments one has coupling constants that describe interactions in the high energy physics model. Instead of evolving the system of equations in time as in the equations of epidemic modeling in the renomalization group equation one has evolution with respect to energy scale whose logarithm is denoted as t. Perhaps the most important system of renormalization group equations of high energy physics is the equations governing the coupling constants of the standard model of high energy physics [9] . These can be written as:

with the right hand side give by the beta functions of the theory whose coefficients are computed using Feynman diagrams and are given by:

with initial conditions given by:

In these equations λ is the coupling constant of the Higgs boson, y t is the coupling constant of the top quark and (g 1 , g 2 , g 3 ) are the coupling constants the electroweak and strong interactions. Looking at the structure of the system of differential equations describing the renormalization group one can see that it describes a system of similar complexity to a five compartment epidemic model. In the following sections we will consider models with eight and twelve compartments of even greater complexity. One can study a simpler form of renoramlization group equations associated with the Higgs boson coupled to the top quark in a model referred to as the Yukawa model [10] . For this model one has two coupling constants and the renormalization group equations are given by:

Now defining as in [10] : 

we have:

with β a λ = 2(12a 2 λ + 2a y a λ − a 2 y ) β ay = 5a 2 y (3.8)

In this form we see that the renormalization group equations of the Yukawa theory have a strong resemblance to the system of differential equations used to describe the SIR model. We may then expect that techniques and visualizations used in the renormalization group equations should be applicable to the equations used in epidemic modeling. Below we illustrate this using the example of flow diagrams and mixing terms in the equations. 

The renormalization group flow diagrams illustrate fixed points in the renormalization flow where the coupling constants do not change in scale. Thus they become fixed at whatever value make the beta functions vanish. In the pandemic context a similar phenomena happens when we reach a point in the pandemic when the number of susceptible, infected and recovered individuals are not changing in time. This can occur if the number of susceptible individuals dropped to a level that the pandemic is extinguished or it can also occur if the number of infected individuals reaches a constant level in the population and then the pandemic virus becomes endemic in the population. In such cases the fixed points of the flow diagram represent points of endemic equilibrium in (S, I, R) space [11] [12] . In figure 2 we show the renormalization group flow for the Yukawa theory considered in [10] . It is clear from the diagram that the flow of coupling constant when the scale changes goes to a fixed point in coupling constant space where the flow is stationary and the beta functions or the right hand side of the renormalization group equations vanishes. In figure  3 we show a similar flow diagram for the SIR epidemic model. In this case we also see that the flow goes to a point where the compartments I and S become stationary when the epidemic has ended after the number of susceptible individuals drops to a level so that the epidemic can no longer be sustained. It will also be important to study flow visualizations over longer time scales where the virus can be come endemic and then recur at a later date.

The contact matrix is a symmetric matrix that indicates the number of contacts between different groups. If the matrix is diagonal then the different groups do not mix. For example if two countries are island countries and have no boat or air travel between them then contact matrix will be diagonal. For the renormalization group one has an analogous concept of operator mixing. This occurs because loop graphs induce new operators that were not present in in initial renormalization group evolution. The matrix of coefficients that describe this operator mixing is the anomalous dimension matrix and is analogous to the contact matrix.

For the renormalization group equation the effect of operator mixing is described by the anomalous dimension matrix γ i,j through the equation:

where C i is the coefficient of an operator O i . These matrices have an important physical effect as operators which have vanishing coefficient at one scale but which develop a nonzero value at a different scale through the mixing with other operators in the theory. For the SIR model the contact matrix plays a similar role. In this case a subset population may have a vanishing number of infections at one time but develop infections through the contact with another segment of the population at a later time. This contact matrix is proportional to the number of contacts between the different populations and is represented in the differential equations of the SIRD model in a similar manner to the anomalous dimension matrix in the Renormalization Group Equations. For the SIRD model the effect contact matrix c ij is given by [13] :

In these equations β i is the contact coefficient between susceptible and infected individuals in group i, γ i is the infection rate in group i and µ i is the death rate in group i and c ij is the contact matrix between different groups. In the SIRD model S i is the number of susceptible individuals in group i, I i is the number of infected individulas in group i, R i is the number of recovered individuals in group i and D i is the number of individuals who died in group i. The SIRD model has been applied to pandemics including the ongoing COVID-19 pandemic in [14] [15] [16] . We will give two examples of the use of the SIRD model in grouping by age and also by country in the next two sections.

The contact matrix between different age groups can be important especially if one age group is more susceptible or has a higher mortality rate than the other. For example the older population may be more inclined to severe outcomes if they have weakened immune systems or secondary health conditions. The other way around is possible too. Older segments of the population may have been exposed earlier in life to a similar pathogen and still retain some immunity while the younger population would not been alive during this time and would be completely susceptible. Still another scenario is that the younger population was vaccinated against a similar pathogen as part of routine child immunizations wheres the older population either never received them or the effect of the vaccines have waned leaving them vulnerable. The SIRD model among two age groups is given by an eight compartment model [17] :

The contact matrix between the two age groups is taken to be:

Other parameters we use are: Where group 1 is the age group 0-49 with a low mortality rate and group 2 is the age group 50+ with a substantially higher mortality rate. The initial population of group 1 is taken to be 60,000 and for group 2 the initial population taken to be 40,000 so that the total initial population of the two groups is 100,000. The initial number of infected individuals in group 1 was taken to be 100 and in group 2 it was taken to be 0. Solving these equations with the parameters and initial condition we see that as the infections in group 1 go down but once contact with group 2 is initiated the infections in group 2 increase and one obtains a second peak of infections in group 2 with contributes substantially to the total number of deaths due to the high mortality in group 2. Often during pandemics segments of the population are kept from the susceptible part of the population through stay at home orders. These segments may be added later to susceptible pool at a later stage in the pandemic when restrictions are eased. This can be modeled with nine groups from 0 to 80+ in ten year intervals as in figure 4 or as two age groups one group under 50 and the other over 50 as in figure 5 . In either case the population with a contact matrix which is diagonal initially and is non-diagonal at a later stage. The results from the SIRD model for the the two age group case is shown in figure  6 and figure 7 . We clearly see the second peak of infection in the older age group due to the the effects of the contact matrix and initial infection in the younger age group. Because of the higher mortality rate in the older age group the model predicts that most of the deaths occur there. For the parameters that we chose in the SIRD simulation with a susceptible population of 60,000 in the younger age group and 40,000 in the older age group and we found that the number of deaths were 86 in the younger age group and 1283 in the older age group for at total of 1369 deaths out of the 100,000 population. This is a dramatic illustration of the effect of the contact matrix between different age populations in pandemic modeling.

It is important to note in the early phase of a pandemic such as from COVID-19 the model parameters are not known accurately. These uncertainties can lead to a range of predictions that depend on these uncertainties. Although there are several parameters in the models we simulate we found the parameter β that indicated the interaction between the susceptible and infected populations to be very important. For the simulation results shown figure 6 increasing β by 10 percent lead to a peak in infection in the younger group increasing from 12,000 to 13,800 and and an increase the peak in the older group from 8000 to 9500. Decreasing β by ten percent lead to a peak in infection in the younger group with a decrease from 12,000 to 9500 and and an decrease in the peak in the older group from 8000 to 7060. For the simulation results in figure 7 increasing β by 10 percent lead to deaths in the younger group increasing from 86 to 89 and and an increase in deaths in the older group from 1283 to 1335. Decreasing β by ten percent lead to deaths in the younger to decrease from 86 to 82 and and an increase in deaths in the older group from 1283 to 1212. Figure 6 : Infections as a function of time for two age groups (0 − 49) and 50+ where there is a nondiagonal contact matrix in the SIRD model describing the interaction between the two age groups during the pandemic. Note the presence of two peaks corresponding to peak infections in both age groups. Figure 7 : Deaths as a function of time for two age groups (0 − 49) and 50+ where there is a nondiagonal contact matrix in the SIRD model describing the interaction between the two age groups during the pandemic. Note the older age group leads to far more deaths than the younger age group due to the higher mortality rate in the older age group.

In this section we consider a simple three compartment model involving three countries involved in virus spreading. In country 1 the virus starts an exponential growth pattern. Travel between country 1 and country 2 takes place, followed at a later date by travel from country 2 to 3. The three country model can be described by the twelve compartment SIRD model:

with contact matrix between the countries given by the three by three matrix which we take to be:

Other parameters we use in the simulation are:

With these parameters and also multiplying the contact matrix with a smoothed out step function to model time delay in establishing contact between the countries we obtain the results for infections in figure 8 and deaths figure 9 . Note that when the first country has contact with the second the infections in country 1 have dropped substantially nevertheless there were enough cases to continue the pandemic into country 2. The same scenario then repeats itself after country 2's infections go down and contact starts between country 2 and country 3. The same step like effect can be seen in deaths in figure 9 where deaths in country 2 pick up after contact is initiated with the initial infection in country 1 with the same scenario continuing between country 2 and 3. Note the three separate peaks in the number of infections which arise when travel restrictions are eased between countries and the stair like structure in the number of deaths as the pandemic moves between the countries. The number of deaths occurring assuming each country has a population of 100,000 is 1356 per country for a total 4068 deaths for the total population of the three countries of 300,000. This is also a dramatic illustration of the effect of contact matrix in this case between countries in pandemic modeling. Again it is important to note in the early phase of a pandemic such from COVID-19 the model parameters are not known accurately. For the simulation results shown figure 8 increasing β by 10 percent lead to a peak in infection in the first country increasing from 19,560 to 23,020, an increase the peak in the second country from 19,630 to 23,030, and an increase the peak in the third country from 19,630 to 23,025 . Decreasing β by ten percent lead to a peak in infection in the first country with a decrease from 19,960 to 18,800 , a decrease in the peak in the second country from 19,630 to 16,000 and a decrease in the peak in the third country from 19,630 to 16,000 as well. For the simulation results in figure  9 increasing β by 10 percent lead to deaths in the first country increasing from 1350 to 1445, an increase in deaths in the second country from 1387 to 1445, and an increase in the deaths in the third country from 1387 to 1417. Decreasing β by ten percent lead to deaths in the first country to decrease from 1350 to 1323, deaths in the second country to decrease from 1387 to 1324, and deaths in the third country to decrease from 1387 to 1324 as well. Figure 8 : Infections as a function of time for three countries whose interactions are described by an nondiagonal contact matrix. Note that although the contacts between country 1 and country 3 are zero country 3 nevertheless becomes part of the pandemic because both countries 1 and 3 interact with country 2. Figure 9 : Deaths as a function of time for three countries whose interactions are described by an nondiagonal contact matrix. Note the stair like structure in the number of deaths as the pandemic moves from country 1 to country 2 and from country 2 to country 3.

There are two possible effects of vaccine waning in pandemic modeling. In the first case a vaccine is developed specifically targeting the virus. Like the seasonal flu vaccine the effectiveness of this vaccine can effectively wane perhaps in short enough time scales so that people will lose immunity to the specific virus during the pandemic and will require a booster to regain immunity [18] [25] . In the second case one sees some efficacy from a vaccine which we call nonspecific to the pandemic virus but was developed for a different pathogen but nevertheless may have some positive effects that result in lower morbidity to the pandemic virus. In references [26] [27] [28] nonspecific vaccines that are being studied including the MMR vaccine, the live attenuated Polio vaccine and the BCG Tuberculosis vaccine. The effects of this nonspecific vaccine could also wane over time. This can also lead to sporadic outbreaks of the original virus it was designed to protect against as for example has been seen in outbreaks of mumps. If the nonspecific vaccine has a positive effect against the pandemic then the effect of vaccine waning may also contribute to age dependence of pandemic morbidity as vaccine waning may reduce the effect of the nonspecific vaccine in older populations who have not received a booster and have a longer time period from inoculation. The waning involved in the nonspecific vaccine could be much longer than the time frame of the pandemic and the main effect would be a contribution to the age distribution of pandemic morbidity.

Here we study a set of two SIRD models to study this effect. In the first model we study a SIRD model with two age groups and no effect of the nonspecific vaccine.

In the second model we consider a SIRD model including a possible effect reducing the morbidity with a greater effect on the younger age group due to possible waning of the nonspecific vaccine. The amount of waning by age group is given in figure 10 and the effect on the death rate with or without the hypothetical nonspecific vaccine that we used in the simulation are shown in figure 11 . The model we used for the second SIRD model is given by:

Then we compare the results of these two models to see the effect of the nonspecific vaccine and the pandemic. We see from figure with the same contact matrix c ij from (4.2). We see that there is an age related effect in pandemic modeling somewhat similar to what is seen in the COVID-19 pandemic. It would be very interesting add more compartments and to represent comorbidities in the different age groups to more closely model the different age populations reaction to the disease. Figure 11 : (Left) Death rate for a population without a nonspecific vaccine for two age groups (0 − 49) and 50+. (Right) Deathrate for a population with a nonspecific vaccine taking into account vaccine waning for two age groups (0−49) and 50+. Note the deathrate is reduced more for the younger age group as the nonspecific vaccine wanes less than the older age group.

We also see from figure 12 that the number of deaths in both age groups decrease with the hypothetical nonspecific vaccine. We chose parameters in the SIRD simulation from (6.3) with a susceptible population of 60,000 in the younger age group and 40,000 in the older age group. We found that the number of deaths decreased from 158 to 86 with the use of the hypothetical nonspecific vaccine in the younger age group. The number of deaths in the older age group decreased from 1351 to 1283 with the use of the hypothetical nonspecific vaccine. Overall the number of deaths in the total population of 100,000 decreased from 1509 to 1369 due to the effect of the hypothetical nonspecific vaccine using the parameters Figure 12 : (Left) Deaths for a population without a hypothetical nonspecific vaccine for two age groups (0 − 49) and 50+. (Right) Deaths for a population with a hypothetical nonspecific vaccine taking into account vaccine waning for two age groups (0−49) and 50+. Note the number of deaths is reduced more for the younger age group as the nonspecific vaccine wanes less than in the older age group. (6.3) of the SIRD model. The effect in the older age group is less because of the vaccine waning varying with the number of years from vaccination. The application of a booster of the nonspecific vaccine to the older age group could be expected to lower the mortality even further as more fatalities are recorded in the older age group where the mortality rate is higher. It will be interesting to examine the projections for pandemic modeling as more data becomes available on the effect of booster for nonspecific vaccines [26] [27] [28] . This type of analysis could also be important for specific vaccines as more data becomes available from clinical trials and measurements are made as to the amount waning that occurs for vaccines designed specifically for the pandemic pathogen [29] .

Finally The number of infections per day remain the same with or without the nonspecific vaccines in the SIRD model we considered. This is because the nonspecific vaccines reduce the severity of the illness and lower mortality but do not prevent infection with the pandemic pathogen as would happen with a specific vaccine targeted to the pandemic virus [29] .

For early stages of pandemics model parameters are not known accurately. To gauge the effects of parameter uncertainties in figure 12 we vary the important parameter β to see its effects on the simulation. For figure 12 (left) increasing β by 10 percent increased the deaths in the simulation without the nonspecific vaccine from 158 to 164 for the younger age group and from 1351 to 1406 for the older age group. Decreasing β by 10 percent decreased the deaths in the simulation without the nonspecific vaccine in the younger group from 158 to 150 and from 1351 to 1277 in the older age group. For figure 12 (right) increasing β by 10 percent increased the deaths in the simulation with the nonspecific vaccine from 86 to 89 for the younger age group and from 1283 to 1335 for the older age group. Decreasing β by 10 percent decreased the deaths in the simulation without the nonspecific vaccine in the younger group from 86 to 82 and from 1283 to 1212 in the older age group.

In this paper we have examined some similar concepts in the renormalization group equations of high energy physics and pandemic modeling based on differential equations. This include the nonlinear behavior of the differential equations, the dependence on the initial conditions and renormlalization group evolution, the similar nature of the contact matrix between compartment groups in epidemic modeling and the notion of operator mixing and fixed points in the renormalization group equations. We put these features together in models with mixing between different age groups, mixing with different countries and the effect of waning of nonspecific vaccines that were designed for different diseases but can lessen the severity of infection and mortality for a pathogen associated with a pandemic. In the future it would be interesting to extend these case studies using more data with respect to different regions, age demographics and age related mortality. It will also be interesting to extend the contact matrix study with different countries as travel is resumed between countries with lower rates of infection. Finally we would like to update the model associated with waning of nonspecific vaccines as well as specific vaccines as more data becomes available from clinical trials about the effect of nonspecific and specific vaccines on the rate of mortality for the pandemic.

Nonperturbative Renormalization Group for the Diffusive Epidemic Process

Renormalization Group Approach to Pandemics: The COVID-19 Case

A Contribution to the Mathematical Theory of Epidemics

Simple models for containment of a pandemic

A model for influenza with vaccination and antiviral treatment

The SIR model and the foundation of public health

Influenza models with Wolfram Mathematica

Self-consistence of the Standard Model via the renormalization group analysis

Renormalization-Group Flows and Fixed Points in Yukawa Theories

Fixed point analysis of the Kermack McKendrick SIR Model

A fixed point theorem for a general epidemic model

Social contacts and mixing patterns relevant to the spread of infectious diseases

Estimating and simulating a SIRD model of COVID-19 for many countries, states and cities

Studying the progress of COVID-19 outbreak in India using the SIRD model

Analyzing COVID-19 data using SIRD Models

Stability analysis of an age structured SEIRS model with time delay

Waning immunity and mocrobial vaccines -Workshop of the National Institute of Allergy and Infectious Diseases

Implications of vaccination and waning immunity

Vaccine waning and mumps re-emergence in the United States

How long do vaccines last? The surprising answer may help protect people longer

Modelling the strategies for age specific vaccination scheduling during influenza pandemic outbreaks

SIR-based mathematical modeling of infectious diseases with vaccination and waning immunity

Modeling the effects of vaccination and treatment on pandemic influenza

Impact of A Waning Vaccine and Altered Behavior on the Spread of Influenza

Could an unrelated live attenuated vaccine serve as a preventive measure to dampen septic inflammation associated with COVID-19 infection?

Can existing live vaccines prevent COVID-19?

BCG vaccine protection from severe coronavirus disease COVID-19

Safety and immunogenicity of the ChAdOx1 nCoV-19 vaccine against SARS-CoV-2: a preliminary report of a phase 1/2, single-blind, randomised controlled trial

This manuscript has been authored by employees of Brookhaven Science Associates, LLC under Contract No. de-sc0012704 with the U.S. Department of Energy.