Comparison of the Evolution of the COVID-19 Disease between Romania and Italy

Abstract

After the outbreak of COVID-19 in Italy, thousands of Romanian citizens who worked in Northern Italy, Spain or Germany returned to Romania. Based on the time-dependent susceptible–infected–recovered—SIR model, this paper compares the evolution of the COVID-19 disease between Romania and Italy, assuming that the parameter value of R₀ in the time-dependent SIR model decreases to R₁ < R₀ after publicly announced restrictions by the government, and increases to a value of R₂ < R₁ when the restrictions are lifted. Among other things, we answer the questions about the date and extent of the second peak in Italy and Romania with respect to different values of R₂ and the duration of the restrictions.

1. Introduction

In December 2019, several cases of pneumonia were detected in the Chinese province of Wuhan. Many common pneumonia and respiratory tract infections were excluded from the investigations, and a new coronavirus infection was discovered and confirmed. The new coronavirus was called SARS-CoV-2, and the disease caused by the virus was COVID-19. The virus is highly contagious and could be very dangerous, especially for the elderly and people with chronic diseases. To reduce the number of infected people countries have approached the problem in very different ways, mostly taking relatively strict measures to prevent the spread, some countries have started implementing measures relatively late, while others had a completely different approach of solving the problem, namely considering heard immunization, without setting restrictions.

At the beginning of 2020, information channels appeared more and more frequently in Romania about a new highly contagious virus that wreaked havoc shortly after its appearance in China. Given the fact that China is geographically located at a considerable distance from Romania, the majority of the population at that time did not consider that it could reach the territory of their country. Unfortunately, the global mobility of the population in the world yield a fast spreading of the virus and the first infections in Romania occurred within just one month. It was soon declared a pandemic by World Health Organization (WHO) officials, and the entire academic community mobilized to stop the spread of the disease. Several factors influence the spread of COVID-19. Predicting the course of the spread of the virus is very important for dealing with the disease and its consequences. To predict endemic and pandemic processes, researchers have developed various models.

Zhang and all [1] used the Poisson model to analyze the evolution of the disease in six Western European countries, and the statistical analysis allowed them to make predictions about the duration, maximum point and total percentage of the population that will be infected. Fanelli and Piazza [2] analyzed the dynamics of the new coronavirus in China, Italy and France. The authors simulated the effects of drastic isolation measures taken in Italy to reduce the infection rate.

Qianying et al. [3] proposed a conceptual model for Covid-19 outbreak in Wuhan, Kucharski et al. [4] estimated how the transmission of Covid19 had varied overtime during January 2020 and February 2020 in Wuhan. Rafiq et al. [5] applied a deterministic model to forecast the spread of COVID-19 in India and Okuonghae and Omame [6] analyzed the impact of non-pharmaceutical control measures on the population dynamics of the new Covid-19 in Lagos, Nigeria using a mathematical model.

Since there are many links between Romania and Italy, we try to analyze in this paper whether the evolution of the virus spread in Romania was similar to that in Italy [7]. The occurrence of the COVID-19 disease in Romania was not due to tourism or direct trade, cultural and economic exchange with China, but to the mostly economic and work related mobility of an important part of the Romanian population throughout Europe, especially to Italy, Spain, Germany or Great Britain. We chose a comparative analysis of the evolution of the spread in Italy and Romania, as Italy was the first European country to fight the major outbreak of Covid-19, it is notorious that many Romanian citizens work, live and travel to Italy, and a considerable amount of statistical data was available for the analysis.

According to the Italian National Institute statistics website [8], many Romanian citizens worked in the areas Emilia-Romagna, Lombardy Trentino Alto-Adige, Friuli Venezia Giulia, that were, in April 2020, heavily affected by the disease COVID-19.

The fact that the first documented case of Romanian citizen infected with COVID-19 was a man of Italian origin travelling to Romania and the next two cases were Romanian citizens traveling to areas in Italy affected by the COVID-19 disease comes in support of this theory of causality

The large number of isolated persons, with whom the statistical data start at 5600 [9,10], is due to the Romanian citizens who entered the country from the time of the first travel restrictions in Lombardy. The number of infected persons continued to increase, either because of the lack of responsibility of the infected persons who did not tell their doctors or friends about their living in the dangerous area of Italy, or because they did not know that they had been infected with the virus. The incubation period (i.e., the time between exposure to the virus and the onset of symptoms) for COVID-19 infected people is estimated to be on average five to six days, but it can be up to fourteen days [11]. The fact that many people, who came into contact with the virus showed no specific symptoms also facilitated the spread of the virus. Consequently, many people were hospitalized and a large number of medical staff was infected, too.

Government measures in Romania were taken as soon as the first cases of the disease were confirmed. The first case was registered on 27 February 2020 and two more the next day. On the 1st March 2020 the terms “quarantine” and “isolated person” were defined. The first data show that the number of people quarantined or in isolation was about 5600 [9,10], after medical check points were set in place at the border for citizens coming from the restricted areas of Italy.

The next restrictions were: a ban on visiting family members in hospitals, a ban on public or private events in open or enclosed spaces with more than 1000 people (8 March 2020), the suspension of the pre-school education process and the suspension of car, air and rail transport with Italy (10 March 2020).

When the number of 200 infected was exceeded (16 March 2020), a state of emergency was declared and other restrictive conditions were imposed: restaurants, hotels, bars would be closed, cultural, artistic, scientific, sports and entertainment events would be suspended, flights with Spain would be suspended and the restriction for travelling to/from Italy extended, universities would go online.

On 21 March 2020, when there were 367 confirmed infected patients, 4207 people in official quarantine and 55198 people in isolation, the restrictions were as follows: the borders were closed, shopping centers were closed, movement of people outside the home was prohibited between 22.00–06.00, dental surgeries, beauty centers were also closed. For private events where liturgical/religious acts are celebrated (weddings, baptisms, funerals), the number of participants had been limited to eight persons.

After statistical analysis of the data from [9,10] using variance and mean value analysis using the tests of Student and Fisher, it was found that the evolution of the COVID-19 infection in Romania differs significantly from that in Italy (p-value 0.001).

The national reports of Romania [12] and Italy [8] for the age distribution of the population reveals that in Romania (Italy) the older population, considered vulnerable to the COVID-19 infection, represents 36.4% (58.8%), respectively.

Italy imposed a restriction on national quarantine on 9 March 2020, when the number of infected persons was 9172, the number of daily infected persons was over 1000 and the number of deaths was 463 per day. In Romania, the maximum restrictions were imposed from 22 March 2020, when the number of infected persons was 433, the number of infected persons per day was 66 and the number of deaths on that day was three. With a population that represents one third of the Italian population, Romania reported 1/15 infected compared to Italy and the number of deaths is less than 1/35 compared to Italy.

Based on real data (see Figure 1) related to the first wave of the COVID -19 pandemic in Italy and Romania, we estimate the peak date, extent and duration of the second wave. For the calculations we use the Mathematica program to interpolate the functions from the time-dependent SIR model. We also study how the duration of the restrictions and the date when the restrictions were introduced in Italy and Romania affect the first wave and how they affect the second wave of the COVID -19 pandemic.

The paper is organized as follows: in the next section we consider the modified SIR model and two interesting theoretical examples. According to the presented model we adjust the parameters to the real-life evolution for Italy and Romania. In Section 3, we present the results of the corresponding SIR models for Italy and Romania. In Section 4, we present some remarks on immunity and finally, we answer in conclusions the questions asked in Section 2.

2. Time-Dependent SIR Model

In this section, the classical SIR model and the time-dependent SIR model are briefly described [13]. In the classical SIR model:

$$\frac{dS}{d\tau }=-\beta IS,\hspace{1em}\frac{dI}{d\tau }=\beta IS-\gamma I,\hspace{1em}\frac{dR}{d\tau }=\gamma I,$$

(1)

where $S\left(\tau \right)$, $I\left(\tau \right)$ and $R\left(\tau \right)$ denotes the number of susceptible persons, infected persons and persons recovered at time $\tau$, respectively, the population of size $N$ is closed, i.e., $N=S\left(\tau \right)+I\left(\tau \right)+R\left(\tau \right)$. The constants $\beta N$ and $\gamma$ stand for: $\beta N$ denotes the (average) number of susceptible persons an infected person comes into contact with (and is possible to spreads the disease), while $\frac{1}{\gamma }$ denotes the average time an infected person takes to recover (i.e., the time of being infectious). Therefore, the quotient $R_0=\frac{\beta N}{\gamma }$ is constant and represents the number of new infections caused by an infected person before he/she recovers. In this paper $R_0$ is considered as time-dependent: $R_0=R_0\left(t\right);$ changing after restrictions are announced.

Rewriting Equation (1) in a non-dimensional form, we define:

$$u=\frac{S}{N},\hspace{1em}v=\frac{I}{N},\hspace{1em}w=\frac{R}{N},\hspace{1em}t=\gamma \tau$$

(2)

to obtain

$$\frac{du}{dt}=-R_{0}uv,\hspace{1em}\frac{dv}{dt}=v\left(R_{0}u-1\right),\hspace{1em}\frac{dw}{dt}=v.$$

(3)

From the dynamical point of view [14,15], the $\left(u,v\right)-$ plane is invariant, $v=0$ is the line of singular points, and $\left(u,v\right)=\left(0,1\right)$ is unstable if $R_0>1$, while all other singular points $v\in \left[0,1\right)$ are stable. Note that in Model (3) all variables are less than (or equal to) one (since representing the percentage of the total population). Note that the exact formula characterizing $R_0$ is $\frac{\beta S_0}{\gamma }$ where $S_0$ is the starting number of susceptible persons. However, if the population $N$ is quite big and there are only few initially infected, the density $S_0/N$ is near 1, thus the approximation $\frac{\beta S_0}{\gamma }\approx \frac{\beta N}{\gamma }$ is possible.

We assume (see Figure 2) the following time dependence for $R_0=R_0\left(t\right)$:

• After introduction of the restrictions/interventions, at time $t_1$, the value of $R_0$ decreases to a certain (lower constant) value $R_1$ (i.e., $R_0>R_1$),
• While the restrictions are in force this value remains constant (for time $t\in [t_1,t_2]$),
• After the end of the restrictions/interventions, after time $t_2$, the value of $R_0$ increases to a certain (higher) value $R_2$ (i.e., $R_2>R_1$),
• The total duration of restrictions is denoted by $\Delta t=t_2-t_1$,
• The transitions $R_0\to R_1$ and $R_1\to R_2$ appear in a time comparable to $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\gamma $}\right.$, which is much less than $\Delta t=t_2-t_1$.

This type of time dependent $R_0$ in a SIR-model causes the evolution of the infected-fraction $v=v\left(t\right)$, as shown in Figure 3. The maximum of the 1st wave is $v_1^{*}$ and appears at time $T_1$. The maximum of the 2nd wave is $v_2^{*}$ and appears at time $T_2$.

By $R_3$ a third value for $R_0$ is denoted (see Figure 4), as one considers again some restrictions to avoid the in-avoidable 2nd peak. Note also that the travels were one of the first things that were blocked by the restrictions, so closed population assumption is legitimate under the regime $R_1$. It is known that even in the restricted period a large number of people was traveling from Romania to Germany, for agriculture. The traveling between countries with the same risk-factor was possible after the end of the 1st wave. Even the traveling from a country with higher risk factor to a country with lower risk factor was possible, but migrations were quite limited at this time and a quarantine was mandatory. We assumed the traveling persons were healthy and tested for COVID-19 and they return in the period close to some (low) multiple of $\frac{1}{\gamma }$. This justifies the assumption of closed population also under the regime $R_2$.

First we adjusted the value of $R_0$ at the beginning from the data [16,17] by simply taking logarithms of $I=I\left(t\right)$ and adjusting the data $(t,\log I)$ for the linear regression line. For Italy we found $R_0\in \left[3.72,5.01\right]$, while for Romania we found $R_0\in [3.56,4.8]$. The lower value corresponds to value $\frac{1}{\gamma }=11$ days, while the higher value corresponds to $\frac{1}{\gamma }=14$ days [11]. Next, we read from the graph $I=I\left(t\right)$ [16,17] the date of the beginning of the restrictions in Italy and Romania (i.e., we read $t_1$). Then we adjust the value of $R_1$ for Italy and Romania so that the maximum value of the first wave $v_1^{*}$ appears at real time $T_1$, which in turn was read from the graph [16,17] according to the real data. For the transition-time of $R_0=R_0\left(t\right)$ from $R_0$ to $R_1$ and from $R_1$ to $R_2$ we assumed a value comparable to $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\gamma $}\right.$, but still in accordance to the actual measures taken by authorities. The value of $R_1$ and $\gamma$ was determined to adjust the evolution of $v=v\left(t\right)$ the to the real data. The values for $R_1$ were checked for values $R_1\in [0.4, 0.9]$, according to [18].

The crucial part of the model is the variation of the value for $R_3$ (see also [19]). Additionally, we vary the value of $t_2$ to get its effect on $T_2$ and on the value of the extent fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$.

In the paper we answer the following questions:

• (Q1) How do the publicly announced interventions/restrictions affect parameter $R_1$?
• (Q2) What is the time delay, $\Delta T=T_2-T_1$, in the evolution of the spread of the disease for both countries?
• (Q3) How does parameter $R_2$ affect the fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ ?
• (Q4) What is the effect of the parameter $R_2$ on the time $T_2$?
• (Q5) How does the parameter $\Delta t=t_2-t_1$ (the duration of the restrictions) affect the fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ ?
• (Q6) How does the parameter $\Delta t$ affect the time $T_2$?
• (Q7) What is the effect of possible re-restrictions on the extent and datum of the second peak?

To answer these questions, we look at the system (3) with the real data output and vary the values of $R_1$, $R_2$, $\Delta t$ for Italy and for Romania. For $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\gamma $}\right.$ we used/assumed the value of 11–14 days, as reported in news. To obtain numerical solutions to Equations (3) we used Mathematica 10.4. In the sequel model (3) with $R_0=R_0\left(t\right)$ as shown in Figure 3.

$$\frac{du}{dt}=-R_0\left(t\right)uv,\hspace{1em}\frac{dv}{dt}=v\left(R_0\left(t\right)u-1\right),\hspace{1em}\frac{dw}{dt}=v.$$

(4)

were adjusted to real data for Italy [16] and Romania [17].

In the time dependent SIR model (4) the corresponding systems of ordinary differential equations are non-linear and therefore it is impossible to obtain the exact solution to a particular initial condition problem. Therefore we use the Mathematica 10.4. command $\mathrm{NDSolve},$ which gives results in terms of interpolating function objects. These interpolating functions are representing the numerical solution to a particular initial condition problem as an approximate function. For the time integration the Runge–Kutta 4th order method was used, while for the interpolation of the solution function the Hermite function was used. We used the following Mathematica code:

$$\begin{array}{cc}\mathrm{In}\left[1\right]:=R0[t\_]:=& \mathrm{Piecewise}\left[\right\{\left\{R0, 0\le t\le N1\right\},\{-(R0-R1/N2(t-N1)+R0,N1<t\hfill \\ & \le N1+N2\},\{R1,N1+N2<t\hfill \\ & \le N1+N2+N3\},\{(R2-R1/N4(t-(N1+N2+N3+N4))+R2,N1+N2\hfill \\ & +N3<t\le N1+N2+N3+N4\},\{R2,N1+N2+N3+N4<t\hfill \\ & \le N1+N2+N3+15N4\left\}\right\}]\hfill \end{array}$$

$$\begin{array}{cc}\mathrm{In}\left[2\right]:=\mathrm{sol}1=& \mathrm{NDSolve}\left[\right\{v^{\prime }\left[t\right]==v\left[t\right]\ast \left(R0\left[t\right]\ast u\left[t\right]-1\right),u^{\prime }\left[t\right]=\hfill \\ & =-R0\left[t\right] u\left[t\right] v\left[t\right],w^{\prime }\left[t\right]==v\left[t\right],v\left[0\right]==r1,u\left[0\right]==1-r1,w\left[0\right]\hfill \\ & ==0\},\left\{v,u,w\right\},\left\{t,0,N1+N2+N3+15 N4\right\}]\hfill \end{array}$$

In the model $r1$ is the initial value of infected fraction $v\left(0\right)=r1$.

Theoretical Example 1. In this example, we show the sensitivity of the output in model (4) on $\Delta t$, with $R2=R0$= 4.0 and R1 = 0.5; i.e., $R_0$ jumps back to its »original« value. The evolution times $t_1$ and $t_2$ under each level are very important. For simplicity, in Figure 4, Figure 5, Figure 6 and Figure 7 we have two slightly different cases of jumps from $R_0$ and the »boomerang effect« looks completely different. In Figure 5, the second peak is larger than the first one, while in Figure 7, the second peak much smaller. We see that the SIR model with the time-dependent $R_0$ is very sensitive to the duration of the restrictions (i.e., to the value of $\Delta t=t_2-t_1$).

Note that the only difference of $R_0=R_0\left(t\right)$ in Figure 5 and Figure 7 are the time intervals (the levels of higher $R_{0H}=4.0$ and lower $R_{0L}=0.5$ is the same).

The time intervals in Figure 4 are as follows: for $t\in$ [0,3.5] we have $R_0=4.0$, for $t\in$ [3.5,5.5] we have the descending line, for $t\in$ [5.5,10] we have $R_1=0.5$, for $t\in$ [10,15] we have the ascending line and for $t\in$ [15,25] we have again $R_2=4.0$.

The time intervals in Figure 6 are as follows: for $t\in$ [0,4] we have $R_0=4.0$, for $t\in$ [4,6] we have the descending line, for $t\in$ [6,10] we have $R_1=0.5$, for $t\in$ [10,15] we have the ascending line and for $t\in$ [15,25] we have again $R_2=4.0$.

• For Italy, the parameters of the model (4) are adjusted to real data are found to be

  $$r1=2.319·{10}^{-7}; N1=1.183; N2=3.630; N3=0.726;$$

  $$N4=1.6;R0=4.28;R1=0.4;\frac{1}{\gamma }=13 days.$$

The parameter $R2$ was varied in the range of [1.2,3]. Parameter value $N2+N3=\Delta t$ was varied in the range of [4.3,6.3]. Note that in the dimensionless model (4) this corresponds to a range from almost two months to more than three months.

For Italy, the time is in units $\tau \times 13$ days, for Romania the time is in units $\tau \times 14$ days.

• For Romania, the parameters of the model (4) adjusted to real data are found to be:

  $$r1=1.5·{10}^{-7}; N1=1.42; N2=2.60; N3=1.33;$$

  $$N4=1.6;R0=4.80;R1=0.83;\frac{1}{\gamma }=14 days.$$

The parameter $R2$ was varied in the range of [1.5,2.5]. Parameter $N2+N3=\Delta t$ was varied in the range of [2,8]. Note that in the dimensionless model (4) this corresponds to a range from one to almost four months.

3. Results

Above best-fit-values with $R_0=R_0\left(t\right)$ presented in Figure 8 and Figure 9 for Italy and Romania are varied slightly (i.e., near the best fit values).

3.1. Italy

The overall evolution of the infected fraction in Italy [20] for COVID-19 is presented in the following Figure 10, Figure 11 and Figure 12.

The estimated time of formation of the 2nd peak is 7.8 months after the beginning (in our model 22 February 2020) of the epidemic. According to [16] our model predicts the 2nd peak about the 15th of October 2020, if $R_2=1.6$. Assuming $R_2=1.6,$ the extent of the 2nd wave is 44.5-times larger than the extent of the first wave:

$$\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.=44.5$$

The 1st peak is at 1.79‰, while the 2nd peak is at 7.99%, if $R_2=1.6$. Figure 13 shows the relation between fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ and $R_2$ in range of $R_2\in \left[1.5,2.1\right].$

Figure 14 shows that $T_2$ tends to 5 months, if $R_2$ exceeds 3 (i.e., the 2nd peak could appear in August). But, if $R_2\in \left[1.6,2.0\right]$ the 2nd peak might appear between September and November.

Varying R2 and N2 in model (4) we generated the results for extent $v_2^{*}$ and time $T_2$. For fitting to linear regression line we used the Mathematica command $\mathrm{LinearModelFit}$ while for fitting to exponential dependence $f\left(x\right)=a+b{e}^{-cx}$ the command $\mathrm{NonlinearModelFit}$ of Mathematica 10.4 was used.

The fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ depends linearly on $R_2$. The (linear) correlation (Pearson) coefficient is found to be $0.999865$. The functorial dependence is found to be $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.=-97.31036+89.0063·R_2$. In

Table 1. Statistics table for coefficients.

	Estimate	Standard Error	t-Statistic	p-Value
1	−97.3136	4.16505	−23.3643	2.67331 × 10−6
x	89.0063	2.06599	43.0817	1.27156 × 10−7

the estimates, standard errors, t-Statistics and the corresponding p-value for coefficients $a$ and $b$ in $y=a+bx$ are presented. The p-value stands for the null-hypothesis that the coefficients $a$, $b$ are equal to zero ($y\leftrightarrow \raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ and $x\leftrightarrow R_2$).

The fitted line and the simulated data are shown in Figure 13.

Next, we fitted $T_2$ in months vs. $R_2$ as follows $T_2\left(R_2\right)=a+b{e}^{-c·R_2}$. We found $a=4.7714696$, $b=242.350437$ and $c=2.80811033$. The fitted line and the simulated data are shown in Figure 14.

We fitted the extent fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ in relation to $\Delta t$ and found $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.=44.8798615-0.09911836·\Delta t$. In

Table 2. Statistics table for coefficients.

	Estimate	Standard Error	t-Statistic	p-Value
1	44.8799	0.0871619	514.903	5.24398 × 10−13
x	−0.0991184	0.015565	−6.36801	0.00141214

the estimates, standard errors, t-Statistics and the corresponding p-value for coefficients $a$ and $b$ in $y=a+bx$ are presented. The p-value stands for the null-hypothesis that the coefficients $a$, $b$ are equal to zero ($y\leftrightarrow \raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ and $x\leftrightarrow R_2$).

The fitted line and the simulated data are shown in Figure 15. The correlation coefficient was found to be $-0.943522$.

Finally, we fitted the date of the second peak vs. $\Delta t$ in terms of $\tau$. The correlation was found to be $T_2=3.6552199083+0.8801832505·\Delta t$. The functorial dependence was found to be $0.999797$ with the statistical table for the estimates, standard errors, t-Statistics and the corresponding p-value for coefficients $a$ and $b$ in $y=a+bx$ as stated in

Table 3. Statistics table for coefficients.

	Estimate	Standard Error	t-Statistic	p-Value
1	3.65522	0.0444303	82.2687	5.02856 × 10−9
x	0.880183	0.00793419	110.935	1.12869 × 10−9

. The p-value stands for the null-hypothesis that the coefficient s $a$, $b$ are equal to zero ($y\leftrightarrow T_2$ and $x\leftrightarrow \Delta t$).

The fitted line and the simulated data are shown in Figure 16.

3.2. Romania

The evolution of the infected fraction in Romania [21] for COVID-19 for $R_2=1.6$ is presented in the Figure 17, Figure 18 and Figure 19. In all figures, $\tau$ represents time multiplied by 14 days; i.e., 15$\tau$ means 210 days or cca. 7 months (after 1 March).

According to [17] and [22] the estimated time of the 2nd peak for $R_2=1.6$ is 7 months after the beginning (in our model 1 March 2020) of the epidemic; i.e., about 1 October 2020. For $R_2=2.5$ the 2nd peak is estimated to about 25 July. If $R_2=1.5$ the extent of the 2nd wave is about $15$-times larger than the extent of the first wave (the 1st peak is at 3.9 ‰, while the 2nd peak is estimated at 6.1%):

$$\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.=15.64$$

Figure 20 shows the relation between fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ and $R_2$ in range of $R_2\in \left[1.5,2.5\right].$

In Figure 21 we see that the 2nd peak could appear also about 1 August 2020, if $R_2\to 3$. Therefore, most likely, the 2nd peak will appear between the beginning of August 2020 and end of October 2020, if $R_2\to 1.3$.

Varying R2 and N2 in model (4) we generated the results for extent $v_2^{*}$ and time $T_2$. For fitting to linear regression line we used the Mathematica command $\mathrm{LinearModelFit}$, while for fitting to exponential dependence $f\left(x\right)=a+b{e}^{-cx}$ the command $\mathrm{NonlinearModelFit}$ of Mathematica 10.4 was used [15].

The fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ depends linearly on $R_2$. The linear (Pearson) correlation coefficient is found to be $0.999273$. The functorial dependence is found to be $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.=-48.095079+42.868254·R_2$. In

Table 4. Statistics table for coefficients.

	Estimate	Standard Error	t-Statistic	p-Value
1	−48.0951	1.27627	−37.6841	2.33102 × 10−8
x	42.8683	0.667915	64.1822	9.61961 × 10−10

the estimates, standard errors, t-Statistics and the corresponding p-Value for coefficients $a$ and $b$ in $y=a+bx$ are presented. The p-Value stands for the null-hypothesis that the coefficients $a$, $b$ are equal to zero ($y\leftrightarrow \raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$ and $x\leftrightarrow R_2$).

The fitted line and the simulated data are shown in Figure 20.

Next, we fitted $T_2$ in months vs. $R_2$ as follows $T_2\left(R_2\right)=a+b{e}^{-c·R_2}$. We found $a=4.409702643$, $b=47.626146393$ and $c=1.852615293$. The fitted curve and the simulated data are shown in Figure 21.

We fitted the extent ratio $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.=a+b·\Delta t$ and found $a=30.490269$ and $b=-0.303047$. In

Table 5. Statistics table for coefficients.

	Estimate	Standard Error	t-Statistic	p-Value
1	30.4903	0.0560693	543.796	2.61013 × 10−15
x	−0.303047	0.011882	−25.5048	2.39391 × 10−7

the estimates, standard errors, t-Statistics and the corresponding p-value for coefficients $a$ and $b$ are presented. The p-value stands for the null-hypothesis that the coefficients $a$, $b$ are equal to zero.

The fitted line and the simulated data are shown in Figure 22. The correlation coefficient was found to be $-0.99542$.

Finally, we fitted the date of the second peak vs. $\Delta t$. The correlation was found $0.99981$. The functorial dependence was found to be $T_2=3.758780098+0.59492099·\Delta t$. In

Table 6. Statistics table for coefficients.

	Estimate	Standard Error	t-Statistic	p-Value
1	3.75878	0.0223477	168.196	2.9797 × 10−12
x	0.594921	0.00473583	125.621	1.7159 × 10−11

the estimates, standard errors, t-Statistics and the corresponding p-value for coefficients $a$ and $b$ in $y=a+bx$ are presented. The p-value stands for the null-hypothesis that the coefficients $a$, $b$ are equal to zero ($y\leftrightarrow T_2$ and $x\leftrightarrow \Delta t$).

The fitted line and the simulated data are shown in Figure 23. The dependence $T_2=T_2\left(\Delta t\right)$ is in units of $\tau$ for both variables.

Theoretical Example 2. In this example, we show the sensitivity of the output in the modified model in case of re-restrictions after the beginning of the 2nd peak, as shown in Figure 24. In Figure 25 and Figure 26 we see that the extent of the 2nd peak is always lower but depends on the time $t_3$ of the reducing of $R_2$ to $R_3$. The evolution for case in Figure 24 and Figure 27 is shown in Figure 25 and Figure 26, respectively.

4. Remarks on Immunity

The coronavirus SARS-CoV-2 has only been circulating in human hosts since December 2019, which means that it is simply impossible to know whether immunity to the disease will last longer than nine months. In the meantime, the results only confirm that COVID-19 patients can maintain the adaptive immunity to SARS-CoV-2 for two weeks post-discharge. Evidence from other coronaviruses (e.g., 2002s SARS 1, or 2012s middle Easter respiratory syndrome) suggests that immunity probably lasts longer than that [23]. The duration of immunity to SARS-CoV-2 by a vaccine (once developed) is estimated to be 6–18 months [24].

However, it is still unclear whether someone with immunity could spread the coronavirus to others while fighting off a second infection. If the immune response were strong enough to crush the virus quickly, the person probably wouldn’t transmit it further. A weaker response that allowed some viral replication might not prevent transmission, especially as people without symptoms are known to transmit the coronavirus [23].

Let us just comment on the evolution in Romania: in the 1st wave, the total number of infected persons was approximately 1% of the population, but, after the 2nd wave the (maximum) total number of infected is estimated to be about 60%. Assuming immunity for approximately 12 months this yields the reduction of the number of susceptible persons to $\frac{4}{10}N$. Possibly the extent of the estimated 3rd wave (again assuming constant $R_2=1.6$) will be at most 3% after approximately $20\tau$ (see green curve in Figure 28). It is possible that due to the continuous (and prompt) reactivation of restrictions already the 2nd wave will not be of such extent, as predicted by the proposed model.

Note that the proposed model must be modified after the 2nd wave, because the regime described in Figure 8 and Figure 9 is no longer valid after the 2nd wave. Most likely the value of parameter $R_0$ (with a slight random error) will remain above $R_0=1$. Such a regime is likely to be in operation in Sweden from the beginning (in February 2020) [9,10].

The evolution of COVID-19 in Sweden (i.e., the function of removed persons $R=R\left(t\right)$) is a strictly increasing function (see [25]; real time evolution of the removed persons $R=R\left(t\right)$), while for Italy the same curve runs almost horizontally between May and August (see [26]; real time evolution of the removed persons $R=R\left(t\right)$).

Finally, note that also in Israel [9] the evolution of the COVID-19 disease (until September 2020) is composed of two waves. According to the shape of $R=R\left(t\right)$ (a large horizontal line between 25 April and 15 June 2020 corresponds to the time period between the end of the 1st and the beginning of the 2nd wave) the evolution in Israel also obeys the proposed model (i.e., the variation of $R_0$). The large difference between the extent of the 1st and 2nd wave in Israel is probably due to larger value of $R_2$ compared to Romania and Italy.

5. Conclusions

Since at the time of the outbreak of the COVID epidemic in Europe, many Romanian citizens returned to Romania en masse from Italy, Spain, Great Britain, etc. when the first wave of the epidemic started in Romania, the first wave was relatively more pronounced in Romania than in Italy. However, the number of deaths in Romania during the first wave was far from being as high as in Italy. The proposed model answers the questions asked in the introduction and forecasts about the same extend of the second peak (maximum at about 8%) for both countries (in relative terms).

Note that the proposed model cannot answer the question when (if) will the epidemics end—it just rises up the question on the relation between the time and extent of the 1st and 2nd peak. Neither the prevalence of the disease can be considered by SIR models. To explain how $R_0$ changes, note that in the SIR model the population is assumed to be closed, and $N$ cannot change, neither the rate of recovery $\gamma$. Therefore, the real change appears on the parameter $\beta$. A simple application of the definition of conditioned probability yields $\beta = pc$, where $c$ is a contact parameter (i.e., the probability of contact between a single infected and a single susceptible) and $p$ is a contagion parameter (i.e., the probability that if a contact occurs, then the susceptible becomes infected). However, $p$ cannot change with restrictions, since it depends on the infectivity of the disease. Therefore, a lockdown-type restriction makes $c$ smaller, since it reduces the contacts. The duration of restrictions (i.e., $\Delta t$) and the social distance approach as considered in the model is reflected linearly in terms of the fraction of the maximal extent and in terms of the expected datum of the 2nd peak. On the other hand, the relation $T_2$=$T_2\left(R_2\right)$ seems to be exponential. The model predicts the influence of the beginning and duration of the restrictions to the extent and duration of the 2nd wave.

The number of total infected in the 2nd wave seems to be much higher than the number infected in the 1st wave (note that the 1st peak is measured in per mils while the 2nd peak in per cents). A reasonable question is of course also how many waves will occur. That depends, of course, on the duration of the actual immunity and on whether/when we will have a coronavirus vaccine. Similar to the assessment of the duration of immunity, the answer to this question should be compared with cases from real-life. From history we know of cases where the virus has decimated the population, so caution is more than appropriate in combating COVID-19. One of the most affected diseases in Europe was Spanish flu. In 1918, the Spanish flu infected 500 million people with the H1N1 flu virus in four waves, which was one third of the world’s population at that time (and killed 17-50 million people). If there are many waves, and if the immunity to SARS-CoV-2 is as strong as the immunity to SARS-CoV-1, the next waves will not be as strong as the 2nd wave. However, this is left for the future work.

To conclude we give answers to the questions (Q1)–(Q7):

• Answer to (Q1): In Italy, the final $R_2$=0.4, however the average $R_2=2.017$, with $\frac{1}{\gamma }=13$ days.
• In Romania, the final R2 = $0.83,$ however the average $R_2$ = 2.14, with $\frac{1}{\gamma }=14$ days.
• Answer to (Q2): The time delay $\Delta T$ varies with $R_2$ and $t_2$, as it can be seen in Figure 16 and Figure 23. Varying of $R_2$ makes sense for the future evolution, while varying of $t_2$, once the restrictions were eased, has just a theoretical meaning. For Italy, $\Delta T$ varies from 2.5 months (if $R_2$ = 3) to 7.3 months (if $R_2$ = 1.5).
• Answer to (Q3): From Figure 13 and Figure 20 we see that the relation is linear and ascending. For $R_2>1.2$ the extent fraction is higher than 10. From Figure 25 and Figure 26 we see that after repeating the restrictions the extent of the 2nd peak can be much lower (i.e., one third) than predicted by Figure 2 and Figure 3.
• Answer to (Q4): The relation $T_2=T_2\left(R_2\right)$ seems to obey the exponential law $T_2\left(R_2\right)=a+b{e}^{-c{R}_2}$. The relation is shown in Figure 14 and Figure 21. The exponential dependence is more considerable in case of Romania (see Figure 21). For higher values of $R_2$ the period between the peaks, $\Delta T,$ is lower. Furthermore, the second peak appears faster and has a greater extent.
• Answer to (Q5): The longer the restriction period, the smaller the fraction $\raisebox{1ex}{$v_2^{*}$}\!\left/ \!\raisebox{-1ex}{$v_1^{*}$}\right.$. From Figure 15 and Figure 22 we see that the relation is linear (descending). This means that prolonging the duration of the restrictions leads to decreasing of the extent of the second wave of infection.
• Answer to (Q6): The parameters $\Delta t$ and $T_2$ are linearly related (ascending). For longer restriction periods one obtains the 2nd peak later. Also, the appearance of the second wave is of lower extent.
• Answer to (Q7): The extent of the 2nd wave in case of re-restrictions is always lower and depends on the time of reacting on the start of the 2nd wave. In the middle of September in Italy, a restriction concerning the closure of nightclubs and the obligatory use of face-masks after 18:00 was introduced. This is a concrete example of lowering $R_2$ to $R_3,$ as noted in the Theoretical Example 2.