(WHO) about a cluster of 41 severe pneumonia cases of unknown aetiology in Wuhan, Hubei province. A newly identified β -coronavirus (whole genome sequenced on 5th of January and isolated on 7th of January) was initially named 2019-nCoV (2019-novel CoronaVirus/12 January 2020) and later renamed to SARS-CoV-2 (Severe Acute Respiratory Syndrome CoronaVirus -2/12 February 2020) by the Coronavirus Study Group of the International Committee on Taxonomy of Viruses [1] . The extensive global transmission led to the declaring coronavirus disease of 2019 (COVID-19) as a pandemic with significant mortality and morbidity [2] . By 27 February 2020, the outbreak of coronavirus disease 2019 (COVID-19) caused 13 285 640 confirmed cases and 578 110 deaths globally, more than severe acute respiratory syndrome (SARS) (8273 cases, 775 deaths) and Middle East respiratory syndrome (MERS) (1139 cases, 431 deaths) caused in 2003 and 2013, respectively. COVID-19 has spread to 216 countries internationally. Total fatality rate of COVID-19 is estimated at 3.46% by far based on published data from China CDC. Average incubation period of COVID-19 is around 6.4 days, ranges from 0 to 14 days. The basic reproductive number (ℜ 0 ) of COVID-19 ranges from 2 to 3.5 at the early phase regardless of different prediction models, which is higher than SARS and MERS [3] . The epidemiological and clinical data collected so far suggest that the disease spectrum of COVID-19 may differ from SARS or MERS. Person-to-person transmission of COVID-19 infection led to the isolation of patients that were subsequently administered a variety of treatments. Extensive measures to reduce person-to-person transmission of COVID-19 have been implemented to control the current outbreak and this has led to a significant reduction in the viral spread and a reduction in new cases. In Bulgaria, the disease began with a limited number of cases. The government took instant, timely and decisive measures, which led to a reduction in the spread of the infection. There was a weak first wave with a peak around April 20, when the new cases reached 91 per a day. There was a clear decrease in morbidity over the next month and a half, of which, due to the implemented measures. After restrictions were lifted and borders were opened, cases began to rise rapidly and right now Bulgaria is in "second" wave, with the intensity of morbidity of around of 250 new cases per a day. A very similar picture could be seen in the rest of Europe, where the mitigation of measures and lifting of the travel restrictions had also led to a rapid increase in the number of positive cases and morbidity rates. A major part of the special efforts to predict, prevent and reduce the transmission of the disease, is to develop a mathematic-computation model for forecasting the COVID -19 spread in Bulgaria in Europe and in the world wide.
1.2. Epidemic Dynamics Modeling and Analysis. The modern approach for studying and forecasting the dynamics of pandemic processes is connected to the use of deterministic mathematical models described by differential equations. They track the quantitative changes over time of main population groups representing different stages of the course of the disease or pandemic (virus carriers, hospitalized, recovered, etc.). The equations in the mathematical models describe the consecution and velocity of transitions from one group to another. The exact dependencies in the equations are controlled by coefficients that take into account the real (measurable) characteristics of the infectious disease and the human population in the respective region and are determined by the empirical data. The most commonly used models, including for predicting the current COVID-19 pandemic, are based on the so-called SIR and SEIR models, which are essentially systems of nonlinear ordinary differential equations. Specifically for COVID-19, see, for example, the Chris Murray Model [4] , the Epidemic Calculator [5] , and the works [6, 7] of Neil Ferguson, Imperial College, London. The same model has been used also to analyze the situation in China [8] . For more detailed and complied analysis see a newly published monograph [9] .
The mathematical and computer modeling of the COVID-19 pandemic is a really interesting and important topic in which mathematical modeling and high -performance computational calculations clearly cooperate. This is a modern tool not only for forecasting, but also for studying the mechanisms of disease spread, for assessing the effect of various interventions or strategies for controlling the pandemic.
In SEIR model the host population is represented by the following compartments:
• Susceptible individuals S(t). These are the people who may be infected and can become virus carriers. Usually at the beginning of the pandemic, as in the case on COVID-19, the whole population of the country is susceptible. • Exposed individuals E(t). These are virus carriers individuals in the latent stage, during which they are not virus spreaders. They usually have no symptoms. • Infectious individuals I(t). These are the individuals who are with strong infectivity virus carriers and virus spreaders. These are the people who may be passing the virus on to others in case of contact. It is supposed that at the start there is at least one infectious individual. • Recovered individuals R(t). These are the recovered from the disease individuals (in the same group we calculate also the deceased individuals).
The total population size is supposed to be a constant N = S(t) + E(t) + I(t) + R(t) during the considered pandemic period.
In SEIR model the quantitative changes over time are described by the following Cauchy problem for system of nonlinear ordinary differential equations (under assumptions that S(t), E(t), I(t) and R(t) are differentiable functions):
(1) SEIR model:
Here S 0 > 0, E 0 ≥ 0, I 0 > 0, R 0 ≥ 0 are the initial numbers of individuals in the fourth compartments. The model involves three parameters:
• Infection (transmission) rate β of the biological epidemic, which is a product of the of contact rate r between people in the population and the probability of transmission β 0 upon contact between an Infectious individual and a Susceptible individual;
where T e denoted average latency time (incubation period). • Rate of the Infectious people transform to the Recovered people, i.e. recovery rate γ = 1 T i , where T i denoted the average recovery time.
A fundamental concept in epidemiology is the basic reproduction number ℜ 0 , which measures the speed of spread of an infectious disease. In SEIR model (1) (see [10, 11] ) it is
When ℜ 0 > 1 the speed of spread is positive and SEIR model gives exponential growth of infectious individuals. If ℜ 0 = 1, then the speed of spread is zero (disease threshold). If ℜ 0 < 1 the spread speed is negative and the infection subsides. SEIR model (1) can be considered as generalization of the SIR model, suggested by Kermack and McKendrick [12] . SIR model has been widely applied in dynamic transmission modeling of a directly transmitted pathogen such as influenza. The latent period is ignored and the compartment of Exposed individuals absents in it:
(3) SIR model:
Here the parameters are the infection rate β 1 and the recovery rate γ 1 . The basic reproduction number ℜ 0 in SIR model (3) is defined as the ratio of the infection rate β 1 to recovery rate γ 1 . The SIR and SEIR models with constant coefficients are useful for modeling of rapid spread of disease to a large number of people in a given population within a short period of time (see [13] ). Seasonal epidemics or controlled pandemic by a switching prevention strategy, which may include vaccination and outbreak control measures are described by models with time-depended transmission parameters (see [14, 15, 16] ). In these cases the number of infectious individuals I(t) can have many local maxima (see Figure 1 (b)), instead of one pick (in the cases of constant coefficients, Figure 1 (a)), and they can be controlled by good prevention strategy. As is mentioned in [8] the value of the infection rate β can be changed by various government control measures. For example, educational facilities closing, gathering restrictions, travel restrictions, cancelation of mass gatherings, and similar control measures impact on the values of of the rate r of contact between people in the population. On the other hand, the recovery rate γ can be changed, for example, by methodological improvement on the diagnosis and treatment strategy. Therefore it is reasonable to use time-depended models to describe the COVID-19 pandemic dynamics.
1.3. Spatial Dynamics. The spatial disease dynamics is among the challenging topics to better understanding of future development of COVID-19 and other related pandemics, due to the ever stronger impact of human mobility to infectious disease dynamics on larger geographical scales. For SEIR model, this means generalization of the basic system of ordinary differential equations to systems of either integral equations, partial differential equations or coupled ordinary differential equations on a weighted graph, thus including the spatial diffusion (see, e.g., [17, 18, 19, 20] ). At a certain level of epidemic development, the long range human mobility is responsible for the rapid geographical spread of emergent infectious diseases, thus violating the Brownian motion hypotheses. Although significantly more complex, the emerging fractional diffusion models provide new opportunities for better analysis of such non-local phenomena. The successful development of spatial time-depended inverse SEIR models will require substantially bigger amount and more reliable data to determine the related coefficients, which are by default strongly heterogeneous in space. The data for time-dependent density of individuals is just one of the challenges. For example, in the numerical tests presented in [17] , data for the flux of dollar bills in US is used to parametrize the model of spatial dynamics. The analysis of spatial dynamics of epidemic diseases is a topic with great social and economic impact. The common understanding of this is rapidly increasing during the ongoing COVID-19 pandemic. This is convincingly supported by the following example. The Institute for Health Metrics and Evaluation (IHME) is an independent research center at the University of Washington, financially supported by the Bill & Melinda Gates Foundation [21] . For the next decade, IHME has formulated the following ambitious target: "Scale up our geospatial analyses toward our aspiration of mapping all diseases, risks, and covariates at 5 km by 5 km resolution or better." In the spirit spatial SEIR models, such (or even coarser) resolution can be achieved only via deep synergy of resources in High Performance Computing (HPC) and Artificial Intelligence (AI) in Big Data environment.
The paper is organized as follows. In Section 2 we introduce a time-depended inverse SEIR model for daily adjustment of the infection, latent, and the recovery rates. Additionally, we propose difference schemes for numerically solving the corresponding direct and inverse problems. In Section 3 this theoretical machinery is used for experimental analysis of the pandemic situation in Bulgaria and for a development of a strategy for mid-time-frame forecasting. Numerical results for three different 14-days-periods are documented and analyzed. Conclusions are outlined in Section 4.
In this section we consider SEIR model with time depended parameters: (4) TD-SEIR model:
The parameters in TD-SEIR model (4) are different in different stages of pandemic and they are specifically for each country. We assume that ω(t) = 1/T e (t) ≤ 1 and γ(t) = 1/T i (t) ≤ 1, where T e (t) and T i (t) are average latency time and average recovery time to the moment t.
One can solve the Cauchy problem in TD-SEIR model (4) knowing the initial data and the value of the parameters. In such way functions describing the quantitative changes in considered compartment of population can be found. But this is not the case. Official data sources relating to the COVID-19 pandemic contains information only for number of active cases and number of Recovered individuals R(t). Since many people who have the disease show no symptoms it is reasonable to assume that active cases are infectious and exposed individuals and will denote them by A(t) = E(t) + I(t). The parameters in the TD-SEIR model (4) are unknown and therefore in order to study the COVID-19 pandemic dynamics the first step is to solve
The inverse problem: To find the parameters β (t), ω(t) and γ(t) in the TD-SEIR model (4) using the knowing real data for the functions A(t) = E(t) + I(t) and R(t).
Actually the confirmed data for A(t) and R(t) are not enough for analytically solvability of the inverse problem. As we have mentioned in Remark 1, to solve the inverse problem in SEIR model we need more information, than in SIR model, which unfortunately is not public. With this reason we fix the latency rate ω(t) to be a constant during the pandemic period. The sensitivity analysis conducted in [8] in case of knowing real data for I(t) (instead of A(t)) shows that the change of γ(t) had strong impact on the infected population in China in contrast to ω(t). Then we are able to find β (t) and γ(t) in TD-SEIR model (4).
We will show how the parameters in TD-SEIR model (4) can be calculated using known real data for COVID-19 disease and solving the corresponding inverse problems. Different procedure for calculation of the parameters in SIR/SEIR models and different data fitting procedures like interpolation splines modeling or least square estimation are used in [22, 23, 24] . We propose different data fitting method.
Remark 2. Let us mention that under the term "real data" we use the official data, given by many public sources. Obviously, we use them instead of others which no one knows exactly! In order to solve the formulated inverse problem we introduced the new functions 
Now we are able to solve the inverse problem by the following steps:
(1) We consider period of n ∈ N days of the pandemic. Let us denote by t 1 , t 2 , . . . ,t n the days of this period. (2) The available dataset contains data for the confirmed Active cases A k = A(t k ) and Recovered individuals R k = R(t k ) at the day t k , k = 1, 2, . . . , n. Therefore the number of Susceptible individuals is also known S k = S(t k ) = N − A k − R k . But the numbers of the Exposed individuals E k = E(t k ) and Infectious individuals I k = I(t k ) are unknown (actually, only the sum E k + I k = A k is known). (3) Let us denote by β k = β (t k ), ω k = ω(t k ) and γ k = γ(t k ) the unknown values of the parameters in TD-SEIR model at the day t k , k = 1, 2, . . . , n − 1. It is naturally to assume that β k , ω k and γ k are nonnegative constants.
(4) We rewrite TD-SEIR model (4) as the following difference problem:
It is natural to assume that at the start of the pandemic all Active cases are Infectious individuals A 1 = I 1 . (5) From the first and the last equations in problem DT D − SEIR k (7) we have
(6) From the third equation in problem T DT − SEIR k (7) we obtain (10)
from which follows 
Remark 3. The obtained values of I β ,k and I γ,k in (9) allow us to calculate the value of the basic reproduction number at the day t k if I γ,k > 0 :
Using TD-SEIR model (6) and the obtained values for parameters β = (β 1 , β 2 . . . , β n−1 ), ω = (ω 1 , ω 2 . . . , ω n−1 ) and γ = (γ 1 , γ 2 , . . . , γ n−1 ) we solve numerically recurrent sequence of initial problems for differential equation systems that corresponded to the days t 1 ,t 2 . . .t n of the considered period of the pandemic. More precisely for t ∈ [t k−1 ,t k ], k = 2, 3, . . . , n, we solve the Cauchy problem (14) TD − SEIR k :
at the end of the day t k−1 are used as initial data in the next problem T D − SEIR k , that describes the quantitative changes over the day t k .
In the next section we apply the method described above to the Bulgarian COVID-19 data and approximate values of Active cases, New cases per day and Recovered individuals. 
To calculate the values of the infection rate β (t) and the recovery rate γ(t) during the considered period with DTD-SEIR model (7) we take average incubation period in Europe T e = 4 (see [25] ) and set ω k = 1/T e . Additionally, we assume that the data, provided by the official Bulgarian Unified Information Portal corresponds to all active cases (the union of exposed and infectious individuals) with a small time shift of T e days. Indeed, most of the time the pandemic in Bulgaria had a clustered spreading profile and most of the daily tested individuals were first-hand contacts of already confirmed cases. Therefore, it is hard to judge when exactly they were infected, e.g., were they only exposed or already infectious T e days ago. The situation with the recovered individuals is similar (one is considered "recovered" after a negative test, but the testing itself might have happen several days after the actual recovery). Only for the dead individuals the exact date is known, but their daily percentage is significantly low and can be neglected. Thus, we decided not to introduce time shifts in the numerical experiments and accept the convention that the official data resembles the active cases and the recovered individuals for the day they were announced. The obtained values β k and γ k corresponding to the days of the considered period are shown on the Figure 2 (a) . The calculated values of the basic reproduction number ℜ 0,k (when γ k > 0) are shown on Figure 2 (b) . To verify TD-SEIR model (6) we solve recurrent sequence of initial problems (14) with obtained values β k and γ k (ω k are fixed) for parameters. On the Figure 3 are shown the real data for active cases A k and Recovered individuals R k during the considered period, and the models curves, obtained using the calculated values X k (t k ) +Y k (t k ) (for active cases) and Z k (t k ) (for recovered individuals). We observe very good agreement between the calculated and the real data. All numerical simulations were performed using Matlab developed by MathWorks [26] . We perform three sets of numerical experiments, related to predicting a 14-day-time-frame in the future, based on the available official data for Bulgaria up to the first day of the time frame. Those time frames are: April 20 -May 03, June 29 -July 12, and July 13 -July 26, respectively. The first time frame is related to the peak of the first wave of the virus. The second time frame is related to the last known date. The third time frame is a forecast for the future, thus we cannot compare the numerical results with official data.
To fix the notation, let the first day of the forecast period be denoted by t K . Then, the last day of the forecast period will be t K+13 . We assume that the information, provided in the morning of every day t k+1 is related to parameter values of the day t k , k = 1, 2, . . . , n − 1, e.g., the 79 reported new cases on June 10 are considered as the C k value for June 9. We always solve the Cauchy problem (14) with initial data
K−1 is performed on two levels. The first level uses the six ratios {β k /β k+1 } K−4 K−9 , respectively {γ k /γ k+1 } K−4 K−9 from the week, before the time frame. The same ratios are used for both corresponding parameters from the first and second week of the forecast, e.g.,
and analogously for the other five ratios. Note, that we denote by β k and γ k the predicted values of the parameters, using the proposed methodology, in order to distinguish them from the corresponding official ground-true data β k and γ k , reported by the Bulgarian officials. We have tried several different approaches for optimal fitting of the parameter ratios, namely various convex combinations of the corresponding ratios for a longer time frame in the past (up to four weeks, see documented data in Table 1 ), but the best results were obtained by simple repetition of the ratios from the very last week, as described above. Possible explanations are related to the rapid changes in the government control measures and the continuous growth of licensed laboratories for PCR testing, which gradually increases the overall tested people every week. The second level is to approximate the accumulative week parameters {β week i , β week i+1 }, respectively {γ week i , γ week i+1 } for the weeks: week i -the days between t K and t K+6 , and week i+1 -the days between t K+7 and t K+13 . For the first week, we have the obvious relation:
For the second week, we fit β week i+1 and γ week i+1 , using the available week information from several previous weeks. We used one fitting algorithm for the first time-frame and another fitting algorithm for the remaining two time-frames, due to particularities of the overall situation in Bulgaria (see Fig. 4 ). Further, once we predict β week i+1 and γ week i+1 , together with the ratios
, we can straightforwardly invert (15) 
Finally, having all the necessary values { β k } K+12 K−1 and { γ k } K+12 K−1 , we can solve the Cauchy problem TD-SEIR k along the whole time period of interest. For the first 14-days-time-frame April 20 -May 03, based on the data in Figure 4 we predicted β week i+1 and γ week i+1 via
We used a parallelogram rule for the infection rate β prediction, due to the strict government control measures at that time (facilities and parks closing, social distancing and travel restrictions, etc.). Thus, in general, the weekly infection rate was expected to monotonically decay towards zero. However, this was not the case for two consecutive weeks with linear growth of the parameter (April 13 -April 19 and April 20 -April 26), which was mainly due to the two big Orthodox holidays in Bulgaria: Palm Sunday (April 12) and Easter (April 19). The strict measures quickly inverted the monotonic behavior of β and the weekly values for April 27 -May 3 were already in vicinity of those before the holidays (from April 6 -April 12). We used a simple averaging rule for the recovery rate γ prediction, due to the early stage of the pandemic and the very long recovery period for infected patient. Therefore, practically most of the recovered individuals, counted within this period, were the deceased ones. This number behaved quite stably and did not oscillate much, meaning that it seemed natural to approximate it with a constant function. The experimental comparison between the official data and our forecast is illustrated on the second row of Figure 5 . We observe good agreement between the predicted and the real data. The number of actual cases (which is the biggest value and the most unpredictable one) is well fitted by our prediction. As expected, the longer the prediction period the larger the error, and our first week prediction is overall more accurate than the second week one. However, we managed to almost exactly capture the daily cases numbers on days 8, 9 and 13, which is quite a long-term forecast for the dynamics of the pandemic. The number of recovered individuals (as the smallest one among the three) is fitted the best.
Using TD-SEIR model (4) at the beginning of April 2020 we made forecast for COVID-19 spread among Bulgarian population. Let us mention that on 12.04.2020 we sent to the Bulgarian National Operating Center (NOC) a forecast for an expected peak the disease (of new cases per day) in Bulgaria on 26.04.2020. Later (at the end of April) the Bulgarian NOC reported a peak in the period 20.04 -26.04. 2020.
For the second 14-days-time-frame June 29 -July 12, based on the data in Figure 4 we predicted β week i+1 and γ week i+1 via
Note that, at the starting day t K of the first and second 14-days-time-frames the weekly behavior of the infection rate β is quite similar. Indeed, we have had strong decrease in the value between weeks i − 3 and i − 2, followed by an increase of size at about half the decrease between weeks i − 2 and i − 1. During the first weeks of prediction, the β week behavior again remains similar: β week i ≈ β week i−3 in both cases. The big difference appears in the second week of prediction. It was already mentioned about the strict government control measures during the first time frame. On the contrary, there were very mild government control measures during the second time frame -Shopping centers, parks, and even stadiums were open, highschool students were taking final exams and were allowed to organize prom balls, etc. As a result, the linear increase of β week holds true up to date. Therefore, we used linear fitting for the infection rate β prediction. On the other hand, the percentage of healed patients among the recovered individuals drastically increased, meaning that simple averaging rule for the recovery rate γ prediction was not adequate, any more. We conducted various numerical experiments and concluded, that the weekly behavior of the infection rate five weeks ago gives the most reliable information about the recovery rate during the current week, i.e., β week i−5 γ week i ≈ const. Therefore, we chose to fit this expression by a constant function. As input data, we used the information from the previous two values, since γ week i is again a prediction and not entirely reliable. The experimental comparison between the official data and our forecast for the second time frame is illustrated on the third row of Figure 5 . We observe an agreement between the predicted and the real data during the first week of the time frame, but not so much during the second one. As in the first forecast, the number of predicted active cases is always larger than the officially documented one, but the daily differences between the two plots are significantly larger than before. Apart from modeling obstructions, this behavior could also be explained with the decreased percentage of critical cases and the decreased average age of infected people. As a result, not all infected people test themselves and the official data In order to analyze the accuracy of the proposed prediction methodology, the relative weekly-based errors for the compartments "Active cases" {A k } K+13 K , "Daily cases"
for both the first and the second 14-days-time-frames are documented in Table 2 . Two different norms have been used -the ℓ 2 one and the sup-norm ℓ ∞ , i.e.,
for the active cases and analogously for the other compartments. We observe good agreement between the corresponding error margins of the two norms, which suggests robustness of the proposed model with respect to the norm choice within the {ℓ p } family. Further, in agreement with Fig. 5 , we observe stable error behavior for the Recovered individuals, independent on the time period and the week of the period. For the remaining two compartments, we confirm that the second week prediction (days {t k } K+13 K+7 ) of each time period is less reliable than the first week prediction (days {t k } K+6 K ) and that in the presence of strict government control measures our predictions are more accurate. Finally, we observe that the number of daily cases is the hardest one to predict, which is expected due to the lower values of that number and its day-by-day fluctuations.
For the third 14-days-time-frame July 13 -July 26 we apply the same setup as for the second time frame, as they are consecutive ones and nothing changed much in between. Here, we do not have official data to compare with, since this time period is in the future with respect to July 12. The forecast can be seen on the forth row of Figure 5 .
In this paper, we applied a time-depended inverse SEIR model for daily adjustment of the infection rate β (t k ) and the recovery rate γ(t k ), based on a priori known information about the number of susceptible individuals S(t k ), the number of the active cases A(t k ), and the number of recovered individuals R(t k ). The problem becomes well-posed, if the latent rate ω(t k ) and the number of infectious individuals I(t k ) is also given. For the former, we assume it to be a constant ω = 1/T e , where the average incubation period in Europe T e = 4 is used. For the latter, we assume that I(t 0 ) = A(t 0 ), i.e., all registered active cases at day one of the pandemic were infectious, and further compute I(t k ), based on all the data from day t k−1 .
Next, we developed a difference scheme for solving the discretized direct Cauchy problem that computes S(t k ), A(t k ), I(t k ), and R(t k ), based on β (t k−1 ), γ(t k−1 ), and ω(t k−1 ). Assuming, once again, that ω(t k ) is constant, we developed a strategy for mid-time-term prediction of all compartments of the host population. The strategy takes into account the level of the government control measures, the new cases/tests dependency on the day of the week, and the length of the healing process. Three sets of numerical tests related to predicting a 14-day-time-frame, based on the available official data for Bulgaria up to the first day of the time frame were conducted. It was observed that in all our experiments the number of recovered individuals was well approximated along the whole time period. The numbers of active and daily cases were well approximated during the first week of the time period, but in the case of weak government control measures and lack of social distancing, the predicted values for the second week significantly exceeded the officially reported ones. A potential explanation is that the truth is somewhere in between, as there are infected people with mild or no symptoms, that do not test themselves, thus, do not appear in the governmental statistics. On the other hand, during strict government control measures with compulsory social distancing, the numbers of active and daily cases can be well approximated throughout the second week, as well.