IMMUNO-EPIDEMIOLOGICAL MODEL-BASED PREDICTION OF FURTHER COVID-19 EPIDEMIC OUTBREAKS DUE TO IMMUNITY WANING

. We develop a new data-driven immuno-epidemiological model with distributed infectivity, recovery and death rates determined from the epidemiological, clinical and experimental data. Immunity in the population is taken into account through the time-dependent number of vaccinated people with diﬀerent numbers of doses and through the acquired immunity for recovered individuals. The model is validated with the available data. We show that for the ﬁrst time from the beginning of pandemic COVID-19 some countries reached collective immunity. However, the epidemic continues because of the emergence of new variant BA.2 with a larger immunity escape or disease transmission rate than the previous BA.1 variant. Large epidemic outbreaks can be expected several months later due to immunity waning. These outbreaks can be restrained by an intensive booster vaccination.

accumulated during these two years, and also the experience of mathematical modelling can be used for the development of more precise models which can be used for effective and reliable predictions.
Conventional compartmental epidemiological models are based on two main assumptions. The first one is that the rate of infection transmission from the class of infected individuals I(t) to the class of susceptible individuals S(t) at time t is proportional to their product. This assumption does not take into account the time-dependent infectivity rate during disease progression. This drawback is partially compensated by the introduction of several other compartments, such as exposed, symptomatic, asymptomatic, quarantined, hospitalized. The second assumption concerns the recovery and death rates that are supposed to be proportional to I(t) at time t and fixed for all the individuals. However, the number of recovered R(t) at time t is determined by the number of infected I(t − τ ), where τ is the disease duration. The distinction between I(t) and I(t − τ ) can be very important during the phases of exponential growth or decay. This difference might be possibly compensated by parameter fitting, but this is not precisely clear and it was not verified. If the disease duration is not fixed for all the infected individuals, then a distributed delay should be considered. More complete compartmental models have similar limitations.
In order to overcome these limitations, we propose a model with distributed disease transmission, recovery and death rates. In other terminology, it is a population structured model, distributed delay model, or immunoepidemiological model. Furthermore, we take into account the population immunity composed by the acquired immunity after disease recovery and vaccine-induced immunity, both of them waning with time. We validate the model with the epidemiological, clinical and experimental data and use it to make prediction of further epidemic progression.
One of the important features of epidemic progression is related to collective immunity. In seasonal influenza epidemics, it is reached 2-3 months after the epidemic threshold followed by rapid epidemic decay [12]. In the case of COVID-19 such epidemic progression was not possible because of relatively high mortality rate. Therefore, the measures of social distancing restrained the epidemic level essentially below the limit of collective immunity. This situation changed with the arrival of the Omicron variant characterized by a higher transmission rate in comparison with the previous variants and less severe disease symptoms with lower mortality rate. Therefore, a large disease outbreak was accompanied by relatively weak social restrictions. Due to reduced social behavior related restrictions together with sufficiently high vaccination level, for the first time from the beginning of pandemic, some countries reached the level of collective immunity corresponding to the existing strains. This is a quite remarkable event, but it did not receive much attention in the scientific literature. We will confirm it below with the modelling results and validation with available data.
We face now a paradoxical situation where collective immunity was reached but the epidemic continues because of the emergence of new variants partially escaping the acquired immunity through earlier infection and/or vaccination. Furthermore, the relaxation of social distancing and of other related measures, immunity waning and stagnating vaccination campaign can lead to strong epidemic bursts several months later. We will present these results below on the basis of mathematical modelling. We describe the basic model in Section 2. In Section 3 we explain the method of determination of time varying epidemiological immunological functions from real data. Epidemic progression described by the model and its validation with real data is given in Section 4. We discuss the obtained results and effectiveness of the proposed model in the concluding section.

Structured immuno-epidemiological model
In this section we introduce a new model with distributed transmission, recovery and death rates. It is an immuno-epidemiological model for which the acquired immunity and immunity due to vaccination determine the effective size of the susceptible compartment.

Basic model
We consider an epidemiological model with four compartments, namely susceptible S(t), infected I(t), recovered R(t) and dead D(t) assuming that their sum is constant, Denote by J(t) newly infected individuals at time t determined by the equality As it is usually assumed in epidemiological models, the number of newly infected individuals is determined by the interaction of infected and susceptible individuals present at time t. However, infection transmission rate between infected and susceptible individuals is not a simple product of the two, it is rather given by the following expression: We assume here that infection transmission rate at time t from the individuals J(η) infected at time η depends on the time difference t − η. In the case of respiratory viral infections, it depends on the viral load in the upper respiratory tract. We will specify the function β(t) below. Under some assumptions, equation (2.3) can be approximated by the conventional formula J(t) = βS(t)I(t)/N . Indeed, if β ≡ const and R, D I, then These assumptions are applicable in the beginning of epidemic. However, in general, they do not hold, and the difference between these two definitions of infection transmission rate can be quite essential in the determination of disease progression. Similarly, the number of newly recovered R n (t) and dead individuals D n (t) can be determined by the following expressions: Distributed recovery and death rates r(t) and d(t) show the probability of recovery and death as functions of time. They are determined from the immunological (clinical) data (see Sect. 3). As before, under some assumptions, this model can be reduced to conventional SIR model [13]. Finally, differentiating equality (2.1) and taking into account (2.2), (2.3), we obtain the equations for S(t) and I(t) as follows: Completing them by equations (2.2) and (2.4), we obtain the formulation of a novel immuno-epidemiological model with distributed infectivity, recovery and death rates.

Vaccination and immunity
In order to take into account the influence of vaccination on epidemic progression, we introduce a new variable m(t) corresponding to the immunity level in the population (cf. [6]). Let us begin with the relation of the immunity level with vaccination. We set where V (t) is the number of vaccinated individuals at time t, and V (t) is the rate of vaccination, the function φ(t) describes how immunity changes with time. It is a positive function with φ(0) = 1 (if vaccine is initially fully efficient), otherwise φ(0) > 0, then it increases up to some maximal value and it decreases after that due to immunity waning. The function m(t) adopts values between 0 and 1.
Note 1: If φ(t) ≡ 1, that is, vaccine does not lose its efficacy with time, then m(t) = V (t)/N . This means that immunity level equals the proportion of vaccinated individuals.
Note 2: If V (t) = 0 for t < t 0 and V (t) = N for t > t 0 , that is, the whole population was vaccinated at time In the case of COVID-19, the second dose of vaccination followed the first one after 2-4 weeks. It is convenient to take their effectiveness into account through a single efficacy function φ 12 (t) (Fig. 1), where the beginning of this curve corresponds to the effectiveness due to first dose and the continuation to the efficacy of second dose. To be specific, φ 12 (t) = φ 1 (t) for t ≤ 21 days and φ 12 (t) = φ 2 (t) for t > 21 days, where φ 1 (t) corresponds to the efficacy of first dose and φ 2 (t) to efficacy of second dose of vaccination. Taking into account the booster dose, we can write Here V 3 (t) and φ 3 (t) correspond to the booster dose. Let us note that m(t) can depend on the virus variant through the functions φ j . Next, consider acquired immunity of the recovered individuals. Then the effective immunity at time t becomes as follows: where the function ψ(t) describes how acquired immunity changes in time. Functions φ 12 , φ 3 and ψ will be determined in Section 3 on the basis of immunological data (neutralizing antibodies). We assume that population groups V 12 , V 3 , and R do not overlap in time. This means that the same person can belong to different classes but not at the same time. In particular, this assumption implies that recovered individuals can be vaccinated when their acquired immunity wanes. Though these conditions are not precisely satisfied, this approximation has a reasonable accuracy. Without this approximation, if the classes overlap, the model becomes essentially more complex, and there are no enough data on the immunity dynamics for such overlapping. Under these assumptions, the population immunity m(t) remains between 0 and 1.
Immunity in the population corresponds to the decrease of the number of susceptible individuals. As such, instead of equality (2.1), we have (2.10) Equation (2.6) remains formally the same but with a different S(t) (as defined in (2.10)): we obtain complete model (2.9)-(2.12) with distributed infectivity, recovery and death rates, and population immunity. The functions V 12 (t) and V 3 (t) in the definition of immunity m(t) are supposed to be given, and they will be determined from the epidemiological data.

Basic reproduction number
Assuming that I = D = m = 0 in the beginning of epidemic (possibly for a particular strain or variant), we can write equation (2.5) in the following form: Equating the terms with the first power of (linearization), we obtain (2.14) Setting λ = 0 at the stability boundary, we obtain the dispersion relation, that is, the condition on the parameters which determines the stability boundary: We assume that β(y) > 0 for 0 ≤ y ≤ τ and β(y) = 0 for y > τ , and define the basic reproduction number Then λ in (2.14) is positive (epidemic growth) if and only if R 0 > 1. Moreover, λ increases with the increase of R 0 . For a positive immunity level m = m 0 > 0, Clearly, R 0 decreases with the increase of immunity level.

Determination of the rate functions
The model developed in the previous section includes infection transmission rate β(t), time-dependent recovery and death rates r(t) and d(t), the rates of immunity waning φ 12 (t), φ 3 (t) and ψ(t), and finally the number of vaccinated individuals V 12 (t) and V 3 (t). All these functions can be determined from the available epidemiological, experimental and clinical data.

Infection transmission rate
We assume that the disease transmission rate is proportional to the viral load of the infected individuals in the upper respiratory tract [9]. Therefore, infected individual becomes infectious even before the appearance of symptoms (incubation period), but the rate of disease transmission is varying in time. Next, we suppose that at the end of the incubation period, when symptoms appear, the person is isolated and gradually stops disease transmission. Hence, time-dependent infection transmission rate takes into account exposed, infectious and isolated individuals.
We represent the disease transmission rate in the form β(τ ) = cP (τ ), where τ is the time-since-infection, P (τ ) is the viral load (cf. [14]), and the proportionality coefficient c depends on the interaction rate between infected and susceptible individuals. It depends on the measures of social distancing (behavioural factor) but not on the virus variant. It will be chosen from the comparison with the epidemic growth curves. Then, using the data given in Figure 1(b) in [9], we determine the function P (τ ) for the Omicron variant (Appendix).
Let us note that the area under the curve (AUC) for the viral load, that is, the integral of the function P (τ ) is associated with the infectivity rate of infected individuals. On the other hand, we obtained in the previous section that basic reproduction number R 0 is determined by the integral of the function β(τ ). Hence, the condition of epidemic growth is determined by the AUC up to a constant factor, independent of virus variant.

Determination of functions φ 12 (t), φ 3 (t) and ψ(t)
The data on the vaccine-induce immunity after the first dose and the second dose are available in [11]. They allow us to determine the function φ 12 (t) (Fig. 1a). In fact, the first two doses can be considered as a single vaccination with time-dependent efficacy rate. Vaccine efficacy follows the antibody level and reaches maximum 15-20 weeks after the second dose of vaccination [11]. The data for antibody level are normalized to the maximum implying that the second dose is fully efficient at its maximal efficacy level. The maximal vaccine efficacy can be slightly less than 100% (depending on vaccine type and trial conditions) but close to it.
We assume that the effectiveness of the booster dose is same as of the second dose of vaccination (that is the acquired immunity level after the second dose), and is denoted by φ 3 (t). The function φ 3 (t) is equal to φ 12 (t) having a little bit different growth pattern from day 1 to day 21. The data for ψ(t), the effectiveness of acquired immunity, after recovery from COVID-19, is collected from [3] and is plotted in Figure 1b.
Vaccination in France started on December 27, 2020. This date is chosen as time t = 0 in formula (2.9). We have collected the time series data of vaccination with first dose (V 1 (t)), second dose (V 2 (t)) and booster dose (V 3 (t)) from December 27, 2020 to February 22, 2022. 1 Combining the V 1 (t) and V 2 (t) from the initiation of first dose of vaccine, we have V 12 (t) = V 1 (t).

Estimation of recovery and death rates for Omicron
We now estimate the recovery and death distribution functions r(t) and d(t) for the Omicron variant. The recovery rate r(t) is estimated by a bimodal gamma distribution using the data given in [25]. Probabilities of recovery and death increase with time and reach their maximums after 7 days for recovery and about 20 days for death. For some individuals, the recovery time is quite long. Available data indicate that bimodal gamma distribution is appropriate for the description of distributed recovery and death rates. We consider the following distribution    (Fig. 2a). The recovery function is defined by r(t) = p 0 F 1 (t), where p 0 is the survival probability, p 0 = 0.9975 (for the Omicron variant).
Since the death rate for the Omicron variant is sufficiently low (0.25%), 2 its influence on epidemic progression is weak. Therefore, we can use the previously obtained distribution without specification of the SARS-CoV-2 variant [13] (Fig. 2b). The best fitted parameter values are as follows:  We can now validate the model with distributed recovery and death rates by the comparison with the epidemiological data. We use equations (2.6) relating newly infected J(t) with recovered and dead. Substituting J(t) from the epidemiological data, 3,4 we find the recovery and death rates. Next, we use these rates to determine total recoveries and deaths and compare them with the data on daily and total recoveries and deaths. Figure 3 shows this comparison for France during the Omicron wave from the beginning of December till the end of February. Recoveries and deaths are in agreement with the data when compared separately (see Fig. B.1 in Appendix B). The parameters of gamma distribution are determined here from the only relevant data from India available for the authors [25] but they are validated by the comparison with daily and total recovery data for France. Country specific data can improve the accuracy of model validation and prediction.
Let us note that epidemiological data in this comparison is taken specifically for the Omicron variant. Though it was dominant during this period of time, some previous variants were also present (basically, Delta variant).

Actual epidemiological situation and prediction for France
In this section, we present the results of numerical simulations of epidemic progression in France during the Omicron wave. We have collected the strain-specific data for the period from December 6, 2021 to February 22, 2022. 3,4 After that, another variant BA.2 becomes dominant in France. We will consider it in the next subsection.
For numerical simulations, we use the functions r(t), d(t), φ 12 (t), φ 3 (t), ψ(t) and P (τ ) estimated in the previous section. The functions V 12 (t) and V 3 (t) are taken from the real data of vaccination in France (Fig. 4a). To estimate the disease transmission rate β(τ ) = cP (τ ), we determine the value c = 1.02 × 10 −5 fitting first 25 days starting from December 6, 2021.
Next, we compare the results of simulations with the available data on the total number of infected by the Omicron variant (Fig. 4b) and on the number of new daily cases (Fig. 4c) for the period from January 1 to February 22, 2022. The number of daily cases passes through its maximal value and decreases. This dynamic corresponds to collective immunity. Agreement between the modelling results and the data validates the model. Further epidemic progression (without new variant) is shown with magenta lines in the same figures. Two cases are considered, without additional vaccinations after February 22 and with new booster dose gradually decreasing to zero during 50 days. In both cases there is a long period without epidemic followed by a strong outbreak in September.

New variant and next outbreaks
In case of emergence of a new strain at time t 0 , the level of immunity m(t) is calculated by formula (2.9) for t < t 0 . For t > t 0 , it is modified as follows: Factor κ takes into account that vaccine-induced immunity and acquired immunity after recovery decreases for the new variant, 0 < κ < 1. Comparison with the epidemiological data in France beginning from the end of February when BA.2 variant becomes dominant allows us to determine κ = 0.825. This value is in agreement with the estimate of vaccine efficiency [28]. Let us note that formula (4.1) does not take into account the transient period between two dominating variants when both of them are present. In order to evaluate how this assumption influences the result, we varied the date t 0 of transition between the two variants. The variation of this date in the model leads to some redistribution of the number of infected individuals between the current and the next outbreaks (Fig. B.2 in Appendix B). However, the overall dynamics and the total number of infected individuals after the second outbreak remains the same. So, we suppose that the approximation of instantaneous transition between the variants is acceptable. More detailed description of a relatively short transition period would lead to a more complex model and would complicate the comparison with the available data.
The results of numerical simulations with model (2.9)-(2.12) are shown in Figure 5. The new growth period due to the emergence of the BA.2 variant will continue during about two months followed by a gradual decay. The next epidemic burst is expected in July-August with the total number of infected individuals during this period about 20 million and about 20 thousand deaths (Fig. B.3 in Appendix B). Since we have not found the data for France about the death rate as a function of time post-infection, we used the data on the number of deaths to determine the parameters of gamma distribution. Once these parameters are determined, we use the model to predict the number of death cases during further epidemic progression (Appendix B).
According to [28], BNT162b2 vaccine is 1.3 times more efficient against BA.1 variant than against BA.2 variant. Efficacy of the same vaccine is studied in [10] for the second dose and booster vaccination. First three months after the second vaccination, it is more efficient against BA.2 (51.7% and 46.6%, respectively), but it is more efficient against BA.1 after the booster dose. Since these data do not allow us to make a definite conclusion about relative efficacy of the vaccine against these two variants, we also model the case where vaccine efficacy is the same for both variants but the disease transmission rate is larger for the BA.2 variant. The value of the coefficient c that determines the transmission rate is found by fitting the data during the time interval of 20 days from t 0 (February 22, 2022). In both cases, a large epidemic outbreak is expected in July-August with approximately the same total number of newly infected individuals (Fig. B.3) in Appendix B.
The main cause of this outbreak predicted by the model is immunity waning. Vaccine-induced immunity of the booster dose in the end of 2021, and acquired immunity after the Omicron variant will essentially wane by the summer 2022 resulting in strong epidemic progression. The model does not take into account possible variation in the measures of social behavior related restrictions.
Additional booster dose is necessary to restrain epidemic progression. Vaccination intensity similar to the booster dose in the end of 2021 will allow complete suppression of the epidemic outbreak (curve 4 in Fig. 5). A similar modelling is carried out for Israel (Fig. 6). The results are qualitatively similar with a possibly larger first outbreak in the spring and smaller outbreak in the summer. The same value κ = 0.825 related to the immunity escape of the BA.2 variant as for France seems to give an overestimated growth in March 2022. So, we present also simulations for κ = 0.88. The first predicted peak decreases but the second one remains almost the same. The case with a constant κ and variable disease transmission rate c is shown in Figure B.4 (Appendix B).

Model
In this work we develop an epidemiological model with distributed disease transmission rate, recovery and death rates. Disease transmission rate is assume to be proportional to the viral load determined in the experiments on cell culture and tissue [9]. The proportionality coefficient depends on human behavior related factors (wearing mask, social distancing, etc.), and it is determined by the comparison with the epidemic growth curves. Time-dependent recovery and death rates are determined from the epidemiological data for daily infection [13,25]. The model is validated by the comparison with the data for daily as well as total recovery and death rates. Next, we introduce population immunity in the model (cf. [6]) which allows us to take into account vaccineinduced immunity and acquired immunity on epidemic progression. Immunity waning in both cases is an important factor which can significantly influence the emergence of new epidemic outbreaks.
Clearly, immunity level in the population depends on the vaccination dynamics. We determine the number of vaccinated people (for each dose) from the corresponding data. Another possible approach considered in some modelling works is to consider it as an unknown variable and to introduce an additional equation in the model. However, the vaccination rate in such equation is an unknown time-dependent parameter which should be taken from the data on vaccination dynamics. Therefore, this is unnecessary model complication leading to an additional intermediate compartment in model formulation.
Conventional compartmental epidemiological models, such as SIR model, can be obtained as some particular cases of the model suggested in this work (see [13]). However, this reduction requires rather strong assumptions which are not, in general, satisfied. In particular, we showed in [13] that SIR model underestimates the growth rate and the maximal current number of infected individuals.
The model is derived under certain simplifying assumptions. Some of them are related to the definition of population immunity (Sect. 3). Furthermore, we consider a homogeneous population without age structure and different social groups, and with the same individual immunity. In particular, different individual immunity can differentiate symptomatic and asymptomatic disease progression. There are various other factors discussed in mathematical modelling and epidemiological literature. Our objective in this work is to introduce a relatively simple but effective model. The essential component of the model, namely the level of immunity, is determined from the clinical and epidemiological data. Its validation with the data allows us to expect that it gives appropriate predictions of epidemic progression. More complete models taking into account additional factors indicated above can be considered in further investigations.

Collective immunity
One of the important aspects of epidemic progression is related to collective immunity. When the number of recovered and vaccinated individuals becomes sufficiently large, the epidemic decays because effective basic reproduction number becomes less than 1. This is the usual situation for seasonal influenza. This scenario was impossible for COVID-19 because it would lead to large number of deaths and to collapse of the public health system. The measures of social distancing, and other relevant restrictions, restrained epidemic progression and kept its level essentially below the limit of collective immunity.
Since the Omicron variant is accompanied by mild disease severity and significantly lower mortality rate, the measures of various social restrictions were not reinforced during its emergence in the end of 2021. Together with a high disease transmission rate this situation has lead to a large epidemic outbreak in some countries. In France, there are more than 10 million cases of Omicron SARS-CoV-2 infection, which is comparable with the prevalence for seasonal influenza epidemics with a similar mortality level.
Numerical simulations presented in Section 4 show that collective immunity was reached in January with further epidemic decay due to acquired immunity (including vaccination) and not due to reinforcement of social distancing, as previously during COVID-19 pandemic. However, a new growth of the epidemic started in the end of February, even before the previous wave was completely finished. The BA.2 variant of Omicron determines this new epidemic wave due to further escape from immunity or increased transmission rate in comparison with the earlier variant.

Reliability of model prediction
In general, the main objectives of mathematical modelling consist in the explanation of the mechanisms governing some process and in the prediction of its future progression. The latter is particularly important in the case of large scale epidemics. In spite of enormous modelling efforts during COVID-19 pandemic, reliable prediction time of mathematical models remains limited to a couple of weeks during which ongoing dynamics of the epidemic can be extrapolated.
There are several factors behind this limitation. First of all, the measures of social distancing were imposed by public authorities and they restricted "natural" epidemic progression. Next, and most important, emergence and infectivity of new variants are impossible to predict. Finally, complexity of modern society with multiple social groups, local and global connections, different vaccination strategies and immunity levels, and so on make the detailed description of epidemic progression very complex and hardly possible.
This situation significantly altered in the case of the Omicron variant. The measures of social distancing and other related restrictions remain more or less the same during last several months with some tendency to relaxation but without drastic changes. The heterogeneity of the society becomes less important in the case of weak measures of social distancing because of the effect of mixing (like a mixture of two different gases can be considered as a homogeneous medium). Different population groups are mixed due to social interaction leading to some averaged characteristics. However, this hypothesis is not verified in the present work. It should be verified in further investigations.
The question about the emergence of new variants should be considered in different time scales. In the medium time interval of several months we can build the model and make some predictions on the basis of the existing variants. As such, BA.2 variant determines the ongoing middle size outbreak. This prediction is in agreement with the model in [5]. However, the next outbreak is expected to be larger due to immunity waning of the booster dose at the end of 2021 (not considered in [5]).
This modelling does not take into account possible emergence of new variants during this time period. If a new highly infectious variant appears during spring and summer, it may not still have enough time to overcompete BA.2. If this new variant can become dominant, it will increase the size of the estimated outbreak. From this point of view, we can say that the model gives the estimate from below.
In a longer time scale, after the end of the summer outbreak, we cannot make a scientifically justified prediction about the emergence of new variants. However, epidemic progression shows well established periodicity clearly seen for the whole-world data with new outbreaks emerging about every four months with smaller or larger peaks. Therefore, we can expect another outbreak in the autumn with another virus variant possibly emerging during the summer outbreak.

Conclusions
The model developed in this work takes into account the main epidemic characteristics and immunological properties at the population level with distributed infection transmission, recovery and death rates, population immunity and immunity waning. The model predicts large epidemic outbreaks during summer 2022 due to the fading of acquired immunity of the booster dose massively applied in the end of 2021 and of the acquired immunity of people recovered from the previous wave due to the Omicron variant. New large scale booster vaccination can restrain expected outbreaks assuming that new variants will not appear in the near future.
We interpolated the viral load (green dots in Fig. A.1) for the Omicron variant by shape-preserving piece-wise cubic interpolation in the interval [0, 3] days. We assume that the viral load decays exponentially after that due  [9], blue curve represents its interpolation, dotted red curve shows the decrease of infection transmission due to isolation.
to isolation. The interpolated function for the viral load is as follows (Fig. A.1):

Appendix B. Additional simulations
Here we present the additional simulation results to better illustrate the model validation and relevant predictions. In order not to overload the main text of the paper with the figures, some of them are given here.

Appendix C. Numerical implementation and parameter estimation
The value of parameters was fitted by minimizing the Sum of Squared Errors (SSE) in such a way that the solutions obtained by the model (2.9)-(2.12) approximate the cumulative number of cases. We applied three methods to minimize the SSE function: first, gradient-based method followed by a step of minimization with a gradient-free method, again followed by a third step of gradient-based method. MATLAB nonlinear leastsquare solvers fmincon and patternsearch are used to fit day-wise cumulative number of cases for the time interval of 25 days in the beginning starting from December 6, 2021. Detailed description of this method and its implementation can be found in [18].
System of equations (2.9)-(2.12) is solved using the fourth order Runge-Kutta method. The integrals are approximated by the trapezoidal rule with time increment 10 −4 . Input functions are estimated in the intermediate time steps by suitable interpolations. One simulation run in Matlab takes 25-40 seconds depending on the duration of time interval.