Mathematical Simulations of the Pertussis Epidemic in England

Abstract

The resurgence of pertussis (whooping cough) becomes a serious problem in many countries including the UK. Differentiation of the accumulated monthly numbers of pertussis cases registered in England in 2023 and 2024 revealed two waves of the epidemic before and after October 2023. Identification of parameters of SIR (susceptible-infectious-removed) model allowed calculating the numbers of infectious persons and reproduction rates. The accumulated and daily numbers of cases and the duration of the first wave were predicted. Since the effective reproduction number is very close to its critical value 1.0, the probably of new outbreaks is very high. May be the, increase of percentage of vaccinated people could decrease this probability.

Keywords

• eol>Pertussis (whooping cough)
• pertussis epidemic in England
• mathematical modelling of infection diseases
• SIR model
• reproduction number
• 1

1. Introduction

The resurgence of pertussis (whooping cough) increases the risk of infant fatality and became a serious problem in many countries including the developed ones [1-16]. Many researchers explain the increasing numbers of cases by decreasing the level of vaccinations and propose to maintain the immunization at the level more than 90% worldwide [4]. Some mathematical studies are trying to improve the control of the pertussis outbreaks [2, 5]. In particular, the age structure of the population was taken into account [2, 5].

A modified SIR (susceptible-infectious-removed) model was used in [5] and some results of calculations are presented. Nevertheless, we have not found any solutions of inverse problems for pertussis outbreaks, i.e. identifications of the model parameters or the reproduction numbers with the use of real datasets as it was done for the COVID-19 pandemic dynamics [17-23].

In this study we will use the accumulated numbers of pertussis cases registered in England in 2023 and 2024 (14 months of a recent outbreak, [1]), SIR model and the method of identification of its parameters, proposed in [24] and successfully used in [17-19, 25, 26]. The numbers of infectious person and the reproduction rates will be calculated. We will try to predict the accumulated numbers of cases and averaged daily numbers of new cases and the duration of the recent pertussis outbreak in England.

2. Materials and Methods

We will use the accumulated numbers Vj of laboratory-confirmed pertussis cases in England, [1] (shown in Table 1). Since only monthly data is available (with different numbers of days), we have calculated the accumulated numbers of days tj starting with 1 January 2023.

Changes in social behavior, quarantine measures; appearing new strains, etc. may change an epidemic dynamics. To detect the corresponding new waves, first and second derivatives of the dV dt t=tj d 2V dt2 t=tj » 0.5ççèae Vt jj++11 -- Vtjj + j ÷

V -Vj-1 ÷ö; t j - t j-1 ø »

1 t j+1 - t j çè t j+1 - t j aeç Vj+1 -Vj - j

V -Vj-1 ÷ö;

÷ t j - t j-1 ø j = 2, 3,..,13; j = 2, 3,..,13; The corresponding monthly characteristics can be obtained by substituting tj by m = j 12t j 365

; i = 1, 2,..,14 in formulae ( 1 ) and ( 2 ). Eq. ( 3 ) takes into account different lenghs of the months.

For every epidemic wave i, the generalized SIR model can be applied, which relates numbers of susceptible S(t), infectious I(t) and removed persons R(t) versus time t, [17, 19]: accumulate numbers of cases could be useful [17, 18]. Since the accumulated numbers of cases Vj are given for every month, there is no need to smooth these values (in comparison with the COVID19 pandemic, where the 7-days smoothing was used for very random daily datasets, [17, 18]). To estimate the average daily numbers of cases and the rate of their change the first and second derivatives can be estimated with the use of simple formulae: dS dt dI dt dR dt

= -a iSI , = a iSI - ri I , = ri I . ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 ) t i » r1i .

Ni = S + I + R Summarizing eqs. ( 4 )-( 6 ) yields zero value of the derivative d (S + I + R) / dt. Then the sum: must be constant for every epidemic wave. We consider the value Ni to be an unknown parameter of the SIR model corresponding to the i-th wave, which must be estimated by the results of observations. There is no need to assume that this constant equals the known volume of population and to reduce the problem to a 2-dimensional one. Many researchers (see, e.g., [5]) use this additional unrealistic condition, which means that before the outbreak all people are susceptible. However, many people are protected by their immunity, distance, lockdowns, etc. In the case of pertussis, the rejection of total susceptibility assumption is especially important, since the disease mainly affects children, most of who are vaccinated. Even in the case of respiratory COVID-19 infection, estimates of the initial number of susceptible people (before the outbreak) in China yielded values between 91 and 138 thousand, (0.006% - 0.01% of the population), [17].

The initial conditions for the set of equations ( 4 )–( 6 ) at the beginning of every epidemic wave ti* can be written as follows

I (ti* ) = Ii, R(ti* ) = Ri , S (ti* ) = Ni - Ii - Ri .

The exact solution of the set of non-linear differential equations ( 4 )-( 6 ) can be obtained using the function

V (t) = I (t) + R(t),

Infection and removal rates ( a i and ri ) are supposed to be constant for every epidemic wave, i.e. for the time periods: ti* £ t £ ti*+1, i = 1, 2, 3,.... The inverse values 1/ ri are the estimations of the average time of spreading infection t i (or the generation time [20]) during i-th epidemic wave ( 7 ) ( 8 ) ( 9 ) ( 10 ) ( 11 ) ( 12 ) which corresponds to the number of victims or the cumulative numbers of cases over time t and has the following form [17, 19]

Fi* (V , Ni , Ii , Ri ,n i ) = a i (t - ti*) ,

V dU Fi* = ò

Ri +Ii (Ni -U )[n i ln(Ni -U ) +U - Ri -n i ln(Ni - Ri - Ii )] ;

r n i = a ii .

Thus, for every set of parameters Ni , Ii , Ri ,n i ,a i , ti* and a fixed value of V , integral ( 12 ) can be calculated and a corresponding moment of time can be determined from ( 11 ). The S(t), I(t) and R(t) values can be calculated with the use of the following equations [17, 19]:

I =n i ln S - S + N - R -n i ln(Ni - Ii - Ri ) ; i i

The derivative dV/dt yields the estimate of the average daily number of new cases and with the use of eqs. ( 5 ) and ( 6 ) can be written as follows:

Using ( 4 )-( 6 ), ( 10 ), and ( 14 ) the second derivative can be expressed as follows, [17]: dV dt

= a i SI d 2V dt 2 = a i2SI ( S - I -n i ) ( 14 ) ( 15 ) ( 16 ) ( 17 ) ( 18 ) ( 19 ) ( 20 )

Zero value of the second derivative ( 15 ) corresponds to the maximum of theoretical estimation ( 14 ) of the average daily (or monthly) numbers of new cases. Corresponding values of susceptible persons Ss and moment of time ts can be calculated with the use of formulas available in [17]. At t>ts the second derivative ( 15 ) cannot be positive, since S(t) diminishes monotonously (see ( 4 )). Therefore, a change in the sign of the second derivative from negative to positive reflects the beginning of a new epidemic wave [17, 18].

The effective reproduction number Rt(t) shows the average number of people infected by one person and can be calculated with the use of different approaches, [21-23, 25, 27, 28]. The generalized SIR model also allows estimating the reproduction numbers for each epidemic wave with the use of the formula, [26]:

Rti (t ) =

S (t) = Ni -V (t) = Ni - I (t) - R(t) n i n i n i

Eqs. ( 15 ) and ( 16 ) illustrate that at the moment ts of maximum daily cases ( 14 ) the reproduction number is still higher than the critical value 1.0 (i.e. Rti (ts ) = 1+ I (ts ) /n i ), but diminishes monotonously and approaches for every epidemic wave i, where the number of susceptible persons at infinity Si¥can be found from the non-linear equation, [17, 26]:

Si¥ = (Ni - Ii - Ri )e

Si¥ -Ni -Ri n i To estimate the final day of the i-th epidemic wave, we can use the condition:

Rti (¥) = Si¥

n i I (tif ) =1.

I (ti0 )=1.

which means that at t > tif , less than one person still spreads the infection [17, 26]. Similar condition can be used to calculate the moment ti0 , when a “zero” patient has appeared [26]:

Usually new epidemic waves are connected with new pathogen strains (see e.g. [29, 30]). Then the value ti0 could indicate the time when the new variant (that caused the i-th epidemic wave) started to circulate in a population. In the case of the first epidemic wave, we can use the values I1 =1; R1 =0. Then t* can be treated as the moment of appearance of “zero” patient. Immediately after 1 an epidemic outbreak, N >> V ³ 1, and integral ( 12 ) can be simplified as follows F1* » lnV / N1.

1 Then eq. ( 11 ) yields

I.e., the exponential growth of the accumulated numbers of cases, the daily numbers of cases dV/dt, and the second derivative ( 15 ) is typical for the initial stages of epidemics. Knowing the value of the parameter g 1 , the time of cases duplication can be calculated as follows [17, 26]:

Different procedures for parameter identification can be found in [17-19, 24, 26], but all of them use the accumulated numbers of cases Vj registered during some period of time. Then the corresponding values Fij*(Vj , Ni , Ii , Ri ,n i ) can be calculated with the use of ( 12 ). Due to the linear relationship ( 11 ) the linear regression between Fij*(Vj , Ni , Ii , Ri ,n i ) and corresponding tj values can be used to find the correlation coefficient r, best fitting lines, and parameters a i and ti* , [31]. Optimal values of the parameters have to ensure maximal values of the correlation coefficient or Fisher function: (21) (22) (23) V = eg1(t-t1* ), g =a N .

1 1 1 t d = where n is the number of observations; m=2 is the number of parameters in the regression equation, [31].

In particular, for the first epidemic wave there are only two independent parameters Ni and n i . After finding their optimal values corresponding to maximum of r, the theoretical SIR curves can be calculated with the use of ( 11 ) - ( 13 ). These dependences correspond to best fitting with the random results of observations and can be used for prediction the future epidemic dynamics (in particular, to estimate the numbers of infectious persons and reproduction rates). These approach was proposed in [16] and successfully used for simulations of a mysterious children disease in Ukrainian city Chernivtsi [24] and different waves of the COVID-19 pandemic [17-19, 25, 26].

3. Results and Discussion

The registered accumulated numbers of pertussis cases [1] and calculated values of monthly derivatives are shown in Fig.1 (according to eqs. ( 1 ) –( 3 ), the daily characteristics listed in Table 1 have to be increased according to eq. ( 3 )). The use of logarithmic scale allows detecting the periods of exponential grows. After the outbreak, corresponding values follow the straight lines (according to eq. (21)). Very good coincidences are visible for Vj values (blue) between February and July 2023; for the first derivative (eq. ( 1 ), black) - between March and June 2023; and for the second derivative (eq. ( 2 ), red) - between March and May 2023.

To estimate the value of the parameter g 1 let us use eq. (21) and values of Vj and tj (j=2 and 7) listed in Table 1.

V7 = eg1(t7 -t2 ) V2 ;

Then (22) and (24) yield 43.34 days as the time of cases duplication. The bacterial pertussis infection spreads much slower, than respiratory ones (e.g., COVID-19 cases duplicated in 2.3 - 3.7 days in February-March 2020, [17]).

After July 2023, the pertussis epidemic dynamics demonstrated a stabilization. The monthly numbers of new cases were almost constant in August, September and October (see Table 1 and the black line in Fig.1), but after October 2023 we see the increasing trend again (see Table 1 and the black line in Fig.1).

The second derivative has changed its sign from negative to positive in October 2023 (see Table 1 and the red line in Fig. 1). It means that a new epidemic wave has started and we cannot use all available Vj values to calculate the optimal values of SIR model. Our attempts to use the algorithm presented in previous Section for complete dataset (n=14), were not successful. Nevertheless, with the use of first nine Vj values, the maximum of the correlation coefficient was isolated.

Fig. 2 represents the calculated SIR curves ( 4 )-( 6 ), the first and second derivatives ( 14 ) and ( 15 ), and the reproduction number ( 16 ). The optimal values of parameters (corresponding to the maximal value of correlation coefficient r=0.998545303649538) are:

Very high values of the correlation coefficient and the Fisher function (F=2401, eq. (23)) demonstrate that the linear dependence ( 11 ) is supported by the results of observations at the confidence level higher than 0.001 (corresponding critical value Fc( 7,1 )=29.2, [32]).

The average time of spreading the infection (according to eq. ( 7 )) t 1 »3.16 days. Since this value is not a duration of the illness, but reflects the speed and quality of isolation of invectives, it could be similar for different diseases in one country. In particular, during the first COVID-19 epidemic wave in the UK, t 1 was estimated as 3.03, [17]. The mean UK household generation time was estimated as 3.2 days for the Delta variant and 4.5 days for the Alpha variant [20]. The value of 21 days used for estimations in [5] looks unrealistic (at least for England)

The value t1* demonstrates that the first pertussis patient has started to infect between January 27 and 28, 2023. The smaller value of V1 (in comparison with the exponential trend for V2 –V7 ) probably reflects the fact that outbreak occurred not in the beginning of January 2023. It would be interesting to compare this result with the observations.

Fig. 2 shows rather good coincidence between the theoretical predictions for the accumulated numbers of cases (the blue line, eq. ( 10 )) and results of observations Vj after the prediction (blue crosses). There are more discrepancies between the theoretical estimation of the daily numbers of new cases (eq. ( 14 ), the black curve) and results of calculations listed in Table 1 (black “crosses”). For the second derivative, the results of observations (Table 1, red “triangles”) significantly deviates from the theoretical red curve (eq. ( 15 )). In particular, the negative value of d2V/dt2 (not shown in Figs. 1 and 2) was registered in September 2023 (see Table 1).

The final number of susceptible person for the first epidemic wave S1¥ = 78,448 (eq. ( 18 )) and the final number of the reproduction number is 0.9649 (eq. ( 17 )), while the initial values were 84,205 and 1.0357, respectively. High S1¥ values and reproduction rates, which are close to the critical value 1.0 (see the magenta curve in Fig. 2), show that the probabilities of a new epidemic wave or a new outbreak (after the final moment t1 f (eq. ( 19 )) are very high. The reproduction numbers in September were still supercritical (varied form 1.0320 to 1.0306). Most likely, we are already observing a new wave after September 2023. Observations during next months will demonstrate the influence of the second wave on the accuracy of long-term predictions.

The presented SIR simulation of the first wave demonstrate rather low numbers of infectious persons (with the maximum of 51 around 9-10 May 2024 and less than 1.0 value after the end of August 2025, see the brown curve). The maximum of the average new daily cases (dV/dt, the black line) was expected around 6-7 May 2024. Nevertheless, a severe second wave can make these predictions irrelevant. Unfortunately, we have only 5 data points corresponding to the second wave (j = 10 - 14). The limited number of observations makes SIR simulations of the second wave impossible, but they can be done later or with the use of recent weekly information about the accumulated numbers of cases.

Accumulated numbers of pertussis cases (blue; curve corresponds to V=I+R, eq. ( 10 )); “circles” to the values taken for identification of SIR model parameters; “crosses” represent observations after the prediction). The average daily numbers of new cases: theory – the black curve, eq. ( 14 ); observations - “crosses”, eq. ( 1 ), Table 1. The second derivative d2V/dt2 : theory – the red curve, eq. ( 15 ); observations -“triangles”, eq. ( 2 ), Table 1. The magenta line represents the reproduction number (eq. ( 16 )). The brown line show the numbers of infectious persons I(t) (eq. ( 13 )).

4. Conclusions

Two waves of the pertussis outbreak in England in 2023 and 2024 were revealed with the use of differentiation of the accumulated numbers of cases. SIR model, the algorithm of its parameter identification, and data corresponding to the period January-September 2023 were used to calculate and predict the accumulated and daily numbers of cases, numbers if infectious persons, and effective reproduction number.

The calculated values of the effective reproduction number are very close to its critical value 1.0. This fact and rather high final numbers of the susceptible persons increase the probability of new outbreaks. May be the, increase of percentage of vaccinated people could decrease this probability.

Conflict of Interest

The author declares no conflict of interests

Ethical Approval Statement

The study does not use any experiments with humans or animals. The data sources are available on the Internet.

Acknowledgements

The author is grateful to Robin Thompson, Matt Keeling, Paul Brown, and Oleksii Rodionov for their support and providing very useful information. The study was supported by INI-LMS Solidarity Programme at the University of Warwick, UK.

References

[1] Confirmed cases of pertussis in England by month - GOV.UK (www .gov.uk)

[2] Rohani P , Zhong X , King AA . Contact network structure explains the changing epidemiology of pertussis . Science. 2010 Nov 12 ; 330 ( 6006 ): 982 - 5 . doi: 10 .1126/science.1194134. PMID: 21071671

[3] Elomaa A , He Q , Minh NN , Mertsola J . Pertussis before and after the introduction of acellular pertussis vaccines in Finland . Vaccine. 2009 Sep 4 ; 27 ( 40 ): 5443 - 9 . doi: 10 .1016/j.vaccine. 2009 . 07 .010. Epub 2009 Jul 21. PMID: 19628060

[4] Galazka A. Control of pertussis in the world . World Health Stat Q . 1992 ; 45 ( 2-3 ): 238 - 47 .

[5] Pesco P , Bergero P , Fabricius G , Hozbor D . Modelling the effect of changes in vaccine effectiveness and transmission contact rates on pertussis epidemiology . Epidemics . 2014 Jun; 7 : 13 - 21 . doi: 10 .1016/j.epidem. 2014 . 04 .001.

[6] Vickers D , Ross AG , Mainar-Jaime RC , Neudorf C , Shah S. CMAJ . Whole-cell and acellular pertussis vaccination programs and rates of pertussis among infants and young children . CMAJ. 2006 Nov 7 ; 175 ( 10 ): 1213 - 7 . doi: 10 .1503/cmaj.051637.

[7] Rendi-Wagner P , Kundi M , Mikolasek A , Vécsei A , Frühwirth M , Kollaritsch H . Hospital-based active surveillance of childhood pertussis in Austria from 1996 to 2003: estimates of incidence and vaccine effectiveness of whole-cell and acellular vaccine . Vaccine. 2006 Aug 14 ; 24 ( 33 - 34): 5960 - 5 . doi: 10 .1016/j.vaccine. 2006 . 05 .011

[8] Celentano L. P. , Massari M. , Paramatti D. , Salmaso S. , Tozzi A. E. EUVAC-NET Group , Resurgence of pertussis in Europe . Pediatr. Infect. Dis. J . 24 , 761 ( 2005 ).

[9] Crowcroft N. S. , Stein C. , Duclos P. , Birmingham M. , How best to estimate the global burden of pertussis? Lancet Infect . Dis . 3 , 413 ( 2003 ).

[10] Crowcroft N. S. , Pebody R. G. , Recent developments in pertussis . Lancet 367 , 1926 ( 2006 ).

[11] Tan T. , Trindade E. , Skowronski D. , Epidemiology of pertussis. Pediatr. Infect. Dis. J . 24 ( suppl .), S10 ( 2005 ).

[12] Carlsson R.-M. , Trollfors B. , Control of pertussis-lessons learnt from a 10-year surveillance programme in Sweden . Vaccine 27 , 5709 ( 2009 ).

[13] M. Keeling , P. Rohani , Modelling Infectious Diseases (Princeton Univ. Press, Princeton, NJ, 2008 ).

[14] von König C. H. , Halperin S. , Riffelmann M. , Guiso N. , Pertussis of adults and infants . Lancet Infect. Dis . 2 , 744 ( 2002 ).

[15] Schellekens J. , von König C. H. , Gardner P. , Pertussis sources of infection and routes of transmission in the vaccination era . Pediatr. Infect. Dis. J . 24 ( suppl .), S19 ( 2005 ).

[16] Broutin H. , Viboud C. , Grenfell B. T. , Miller M. A. , Rohani P. , Impact of vaccination and birth rate on the epidemiology of pertussis: A comparative study in 64 countries . Proc. Biol. Sci . 277 , 3239 ( 2010 ).

[17] Nesteruk I. COVID19 pandemic dynamics . Springer Nature , 2021 , DOI: 10.1007/ 978 -981-33- 6416-5

[18] Nesteruk I. Detections and SIR simulations of the COVID-19 pandemic waves in Ukraine . Comput. Math. Biophys . 2021 ; 9 : 46 - 65 . https://doi.org/10.1515/cmb-2020 -0117

[19] Nesteruk I. Simulations and predictions of COVID-19 pandemic with the use of SIR model . Innov Biosyst Bioeng , 2020 , vol. 4 , no. 2 , 110 - 121 , doi:10.20535/ibb. 2020 . 4 .2.204274. http://ibb.kpi.ua/article/view/2042748.

[20] William S Hart , Elizabeth Miller, Nick J Andrews , Pauline Waight, Philip K Maini, Sebastian Funk , Robin N Thompson. Generation time of the alpha and delta SARS-CoV-2 variants: an epidemiological analysis . Lancet. Infectious diseases. Volume 22, ISSUE 5 , P603 -610, May 01, 2022 . https://doi.org/10.1016/S1473- 3099 ( 22 ) 00001 - 9