A Continuous Model of Three Scenarios of the Infection Process with Delayed Immune Response Factors

Abstract

The course of an infection was modeled as a controlled nonlinear process. Understanding the substantial differences observed in the trajectory of the disease caused by the new coronavirus SARS-CoV-2 is of critical importance at the moment. Numerous factors have been considered to explain the fact that symptoms vary highly among different people and the infection transmission rate varies among local populations. Virus replication within the host cell and the development of an immune response to virus antigens in the body are two interdependent processes, which have aftereffects and depend on the preexisting states of the cell and virus populations. Different scenarios with the same input parameters are necessary to consider for modeling the properties of the states. The efficiency of the immune response is the most important factor, including the time it takes to develop defense responses from three levels of the immune system, which is a complex system of the body. A computational description of infection scenarios was proposed on the basis of a delay differential equation with two values of the time lag. In the new model, transitions between phases of infectious disease depend on the initial virus dose and the delayed immune response to infection. A variation in the dose of the virus and response time can lead to a transition from an acute phase of the disease with overt symptoms to a chronic phase or fatal outcome. Asymptomatic transmission of viral infection was calculated and described in the model as a situation where the virus is rapidly and efficiently suppressed after a short replication phase, while still persisting in the body in minor amounts. An analysis of the model behavior is consistent with the theory that the initial number of virions can affect the quality of the immune response. The reasons that high individual differences are observed in the trajectory of COVID-19 disease and the formation of types of the immune response to coronavirus are still poorly understood. Known trajectories of hepatitis C virus (HCV) infection were used as a basis for model scenarios.

INTRODUCTION

This article further develops the previous works on modeling the extreme and special transition modes of existence that an aggressive species takes in a new environment [1]; however, another mathematical apparatus is used here to study this important problem. The method to be employed in modeling should be constructed and chosen as being the most adequate to the specific biological collision. What should remain unchanged is the principle that the logic of situation development is considered as a set of computational scenarios. A single scenario is rarely possible for the processes of interactions between species [2]. Invasion dynamics substantially varies when foreign species enter various habitats [3]. Likewise, the propagation rates of infectious diseases vary among subpopulations. Infection with a pathogenic virus can be considered as a particular case of invasion processes that occurs in conditions of intense counteraction, which determines the outcome of the interplay between the defense systems of the body and the parasite evolving in the body.

The objective of this work was to provide groundwork to a phenomenological model that describes the differences observed in disease trajectory as dependent on the initial virus dose at infection and the time it takes to develop an efficient response. A continuous model is therefore better suited for explaining the outcomes of fighting a virus compared with a model of iteration chaotic dynamics of the trajectory with a set of controlled bifurcations, which was previously used to describe forest pest irruptions. A delayed regulation of the formation of a response is formalized in the equation with regard to the subthreshold viral load at which the effect of infection on the body becomes a significant factor. The subcritical concentration introduced together with two values of the response lag should vary among different organisms to explain the variation in disease scenario among cases with different initial infectious doses. The model was developed to describe the scenarios by which hepatitis C virus (HCV) infection becomes chronic after an acute phase. An asymptomatic disease development with persistence of the causal agent was considered in the model as an important feature of the unusual course of the new COVID-19 disease.

Given the pandemic, it is interesting to discuss the reasons that the infection transmission rate is so variable and the reasons that several variants are observed for the waning phase in local epidemics. The issues are of particular interest in terms of the lockdown efficiency. The first part of the article analyzes the current data on the qualitative dynamics of the SARS-CoV-2 epidemics and provides several explanations of the differences from the viewpoint of the indefinite regulation of processes. Transition types and a mathematical classification of scenarios observed in particular regions are noted as unusual issues that should receive special attention when modeling the period for the epidemic process to end and predicting a repeated spreading of the new RNA virus, possibly a mutant variant.

DIFFERENT AND COMMON FEATURES OF EPIDEMIC PROCESSES UPON THEIR GLOBALIZATION

Epidemics are long known to follow different scenarios in different countries. As an example, substantial differences between African regions were observed for Ebola virus disease [4]. Certain features and unusual phenomena are characteristic of the current pandemic, giving more room for speculation. Once the SARS-CoV-2 pandemic is over, analyses will continue in epidemiology and mathematics to evaluate the lockdown efficiency in various countries, to compare graphs, and to construct models and prognoses for future epidemics. Opposite conclusions will certainly be made, as was observed in many other cases [5].

Compared with the epidemic wave in the United States, the number of cases increased more slowly by one order of magnitude in Mexico. The worst form of a fast-developing epidemic with an extended Π‑shaped peak is observed in the United States based on general data. The qualitative epidemic dynamics varies greatly among individual states, and it is therefore difficult to obtain well-grounded predictions of the disease incidence for neighbor states.

The apparent variation in virus spreading rate is so great that they joke in Mexico they are ready to build the Mexico–United States border wall themselves. The virus proved to be powerless in Vietnam; i.e., not a single lethal outcome has been recorded in Vietnamese statistics by the time of writing this article. After doing 180 000 tests, Vietnam declared April 20, 2020 that the threat had been warded off and lockdown restrictions ended. National behavioral stereotypes should not be ignored when making predictions. It is possible that Italians suffered more severely because they are communicative and talkative [6], and there are physical grounds for this assumption. A study of virus transmission routes showed that swirls of microscopic saliva droplets are emitted during intense conversations and that the louder one speaks the more virions are emitted and the more likely virus transmission is [7].

A classic example is provided by different consequences that the Spanish flu had in two American cities 1 century ago. Two main scenarios are demonstrated by epidemic (mortality per 100 000 population) plots constructed for Philadelphia and St. Louis. A dichotomy observed in the development of the epidemic processes in 1919 provided a basis for mobility restrictions introduced in March. Isolation measures were taken in a timely manner and all public events were canceled in St. Louis, thus preventing an epidemic outbreak with a Λ-shaped mortality peak. Local scenarios that differ so greatly are impractical to combine in statistics, and aggregate data on disease incidence are confusing. After 2 months of the epidemic, a single Λ-shaped peak did not seem that awful in view of economic collapse. The epidemic curve suddenly started to grow again after flattening in many American cities. The coronavirus transmission rate increased again after a deceleration, thus contributing to further differentiation of epidemic scenarios in 2020.

A crucial question is what factors determine the dichotomy of epidemic scenarios [8]. Either rapid growth occurs in the first 2 weeks and the resources of the health care system fail to cope with the number of severely ill people, or the process takes a smooth (or sporadic) form. We note that Philadelphia with its classic Λ-shaped peak is a seaport and a communication center. Therefore, substantial differences in the initial conditions for triggering an epidemic presumably preexisted between the two cities. Each primary spreader induces a chain of transmissions, which are unpredictably wide in scale in some cases.

A greater number of infected people and possible latent carriers (patients zero) naturally arrive to a seaport from centers of infection, and the infection transmission rate may start accelerating faster without leaving time for proper measures to be taken. The plots shown in Fig. 1 are often referred to as an illustrative example, but it must be admitted that the initial conditions were different in the two cities at the starting point. It is not known for certain where the virus emerged in 1918 and whether a cook in a Kansas military camp was the first victim. After peaking again in large cities with a greater mortality rate, the Spanish flu mysteriously disappeared from Europe in June 1919 forever [9], although the Blue Death was expected to return in 1920. Predicting the time for an epidemic to end is an intricate multifactorial problem for deterministic models.

Fig. 1.

Comparison (by mortality rate) of the development of the Spanish flu epidemics in 1918 between two American cities.

In the United States, a critical situation with the coronavirus was observed in New York in April 2020. New York is the largest center of the Chinese diaspora and a historical emigrant gateway of North America. The new virus will spread everywhere sooner or later as a result of globalization. Many virus spreaders may arrive at once in some towns. Transmission chains will multiply, and then local hospitals will face a series of epidemic waves.

A well-grounded theory [10] states that an epidemic (as well as irruption of an invading species) rapidly accelerates once the affected patient concentration has reached a certain threshold in the population and similarly rapidly wanes once sufficient herd immunity has formed in people recovered from the illness and unsusceptible people. The idea is mathematically based, in particular, on modifications of Whittle’s threshold theorem [11] with an epidemics being considered as a branching stochastic process [12] where everyone has a probability of node transition from the ill state to the immunized state [13].

Once the incidence reached 500 cases per day, the epidemics, in fact, dramatically accelerated in the United States as judged from a renewable plot. A fourfold increase in daily new cases was observed instantly. The increase is seen on the plot as a jump over a threshold, where a certain chain of events creates the critical concentration. In fact, the respective people were ill long ago because there is an indefinite lag in how the plot of daily new cases reflects the state of the epidemic process (with incomplete data). An actual threshold state is the infected case concentration at which an eruptive phase is triggered in the infection site; this state has been passed slightly more to the left on the time axis at a lower number of (detected) new infection cases.

As is seen from a renewable plot of current statistics by country (e.g., the one on the Financial Times web page), lockdown measures were efficient when introduced after detecting single cases, as was done in South Korea. When stay-at-home orders are introduced later and are not strict, they fail to quickly change the quality of the process. The epidemic process appears to be an inertial process. Ireland introduced the stay-at-home orders on March 20, 2020 and made them stricter on March 24, 2020. A plateau was reached in Ireland by April 8, 2020. There were 380 daily new cases according to a smooth logistic curve, but then the epidemic rapidly accelerated, there were 500 new cases on April 9, 2020, and a psychologically important level of 1500 daily new cases was surpassed on April 4, 2020. Ireland provides an interesting example where a model equilibrium value, at which the curve is expected to stop growing (a plateau), turns to be a threshold and unexpectedly leads to a dramatic acceleration of the epidemic. By May 5, 2020, daily new cases did not decrease in Ireland to the level observed on April 8. The ultimate mortality per 100 000 population will probably not differ between Ireland and lockdown-free Sweden. The analysis leads to the hypothesis that a peak is unavoidable when lockdown is delayed; i.e., the principle “if you are late, then you are late” holds.

Sweden did not implement stay-at-home measures, and a classic Λ-shaped peak was not observed; however, volatile dynamics with alternating periods of growths or decreases in daily new cases will be observed for long. On April 21, 2020, there were 700 new cases in Sweden and 3500 in Italy, where “stai a casa” orders were implemented long ago. When qualitatively analyzing the epidemic process, disease incidence data on metropolises should be considered separately from aggregated statistics on the total country. This is especially important in the case of Moscow, New York, and Mexico, where the population density is far greater than the density averaged over the respective country. Cases omitted earlier were recorded at a later time in some countries (France), and the delay distorts the scenario of the process and creates false maximums in the infection statistics.

China’s anti-epidemic measures work only in China. Hong Kong and Singapore have already reported a victory over the epidemic, but Singapore’s authorities expect repeated outbursts because of the need for foreign workers. The curve of new daily cases often becomes alternating after peaking; the variation in detection lag acts as a factor. Patients are most commonly identified when they have already developed symptoms and an immune response. Latent development of the infection process in the body plays an important role and creates the effect of an integral delay in transmission, and a temporal decrease in daily new cases should therefore not be interpreted as a success. In Iran, an apparent second wave of the epidemic has been reported from July after a quiet period.

A smooth scenario of the epidemics (the second derivative should be the focus of attention here) was observed in Germany compared with Italy and Spain, where 16 and 20 days elapsed, respectively, from the start to the first maximum. The scenario with a rapid peak is similar to a threshold-dependent triggering of an outburst after a long-term stabilization in the vicinity of an unstable equilibrium [1] (see Fig. 2). A turbulent mode develops after a plateau in this case.

Fig. 2.

Modeling of the threshold scenario of a population outburst by damping (reproduced from [1]).

In the case of an epidemic, a similarity is observed in some cases for the phase of a threshold-dependent transition from fluctuations to a rapid growth. The end of an epidemic is similarly rapid only in rare cases, like a drop in the population size of insects after their intense propagation.

COVID-19 spreads faster than swine flu and has a latent period that averages 5.5 days [14] and ranges from 0 to 15 days (according to Wuhan data [15]). It was observed as early as in Wuhan that virus transmission is possible at the asymptomatic stage, when infection develops and evades adaptive immunity. The time to fever as a main nonspecific response of the body to respiratory infection (and the time to a release of pyrogenic interleukins) in the case of COVID-19 is greater on average [19] than in the case of the related SARS virus of 2003. Neutralizing antibodies were detected using pseudovirus particles in 79, 92, and 98% of samples collected, respectively, from infected healthcare providers who did not require hospitalization 13–20, 21–27, and 28–41 days after the onset of symptoms [17]. The response was certainly slow compared with responses to the majority of the viruses that cause acute respiratory infections. Data on killer T cells are unavailable. Activation of immune system mediators in response to SARS-CoV-2 is delayed in a substantial part of infected patients because of an important, but still unclear factor. This is advantageous for the virus and allows its latent virus transmission. However, the apparent slow development of the response might protect the body from a greater risk when the initial dose is not considered to be critical.

The issue of what biological factors are responsible for the great variation in latent disease phase, which is unusually extended for an acute respiratory infection, will be investigated and discussed in medicine for many years. Attention is paid to the unusual concomitant symptoms, such as loss of smell and loss of taste, which are uncommon for acute respiratory viral infections [18]. It seems reasonable to assume that symptoms become overt when a certain threshold virion concentration is exceeded. Surprising results have been published by Italian healthcare specialists [19]. The majority (90 and 80%) of a municipal population was tested twice, and 40% of the cases were reported to be asymptomatic. The viral load was found to be much the same in virus carriers with overt symptoms and asymptomatic cases. This observation still lacks a reliable explanation. The term “slow,” rather than “asymptomatic,” infection scenario is safer to use. Mathematically, this means that virus replication persists in stable equilibrium (for an indefinite period of time). If symptoms have not developed within a week, this does not necessarily mean that they will not develop after a year. Certain viruses (e.g., HIV) are known to persist in a subclinical state for many years. Other mechanisms that fight infection may come to exhaustion, or the virus may evade the mechanisms as a result of antigenic drift in the body. The variation in disease progression rate in the case of HIV is explained by the fact that a greater number of viral glycoprotein variants bind with antibodies in some people, while avidity for a minor portion of virions is observed in some others. Immunity responds to antigens as individual regions of proteins of an infectious agent, but these regions can be variable. Some regions change rapidly, while others are conserved. Quantitative antibody-related parameters of the immune response do not characterize the efficiency of neutralization of a virus when stably affine immunoglobulin variants are lacking. This is a probabilistic factor in response efficiency. Cases are known where virus-induced systemic inflammation persists after the virus has been eliminated.

In 2020, modeling specialists received unique opportunities to monitor the epidemic trajectories of various countries in real time and to make predictions that can be verified quickly. However, the predictions have failed one after another because of the specific properties of COVID-19. Data on percent positive tests in the total daily test number would be informative in the current situation, but are not available. Long-term model predictions are not feasible to make for the future situation with COVID-19 until the main questions are answered: how immunity forms, how rapidly a response arises in different components of the immune system, and how long recovered patients remain protected. Data on the very early stage of COVID-19 development would be interesting to have for modeling because transition modes and threshold effects may depend on the initial conditions in nonlinear systems. Another important question is whether the initial infectious dose affects the rate of pneumonia development and the immune response, which may be aggressive as a nuclear chain reaction and ineffective at the same time. A sensitivity to N(0) is essential for an adequate model of the variants of new disease development.

FEATURES OF HOW THE STATISTICS OF THE EPIDEMIC ARE PERCEIVED

Available statistical data on virus spreading does not guarantee correct conclusions. Systematic testing of the total population was immediately launched only in Iceland and, in part, Australia, once it was understood that an infected person starts spreading the virus as early as at 3 days of infection. Data from these countries will be used as reference statistics in future analyses of COVID-19. It was early understood on the island that infected, rather than healthy, people are those who need isolation. Funds are more effectively spent on means for mass testing. Figure 3 shows a comparison with the percent distribution of the coronavirus statistic of the Netherlands, were tests were performed only for patients hospitalized with symptoms of viral acute respiratory infection. If this fact is not known, a biologist or statistician may decide in the future that some genetic differences in virus susceptibility have accumulated in the Icelandic population over its long isolation on the island with low genetic diversity.

Fig. 3.

The relative numbers of coronavirus-infected cases by age group in Iceland and the Netherlands (from https://www.covid.is).

Data from Iceland demolish the myth that elderly people are more often affected by the coronavirus. Middle-aged people are actually affected more often. Complications are more intense in the elderly on average. Given this fact, a maximum would indeed be observed at 60+ if tests were performed only in those seeking medical advice. COVID is not the first infection that is milder for people at a young age. Many are aware how unwelcome rubella is in adulthood. Our immune system is a self-regulating taught system that protects the body and consists of several levels. It is easier for the immune system to detect and fight new agents at a younger age than to quickly retrain itself in the face of sudden threats later in life. When a simple explanation is possible, the most logical connections between facts are better for use in mathematical modeling.

PROBLEMS IN PREDICTING AND COMPARING THE SEQUELS OF EPIDEMICS

The basic reproduction number R₀ is estimated for a virus causing an epidemic. The number shows how many other people get infection from one carrier on average. In the case of the SARS-CoV-2 virus, R₀ was estimated at 2.0 by mathematical modeling [20] and 2.5 by analyzing the Wuhan statistics [15] (there are grounds to believe that the data were underreported). The R₀ values are statistics for a reference book. The characteristic fails to explain rapid local outbursts. An apparently higher basic reproduction number was observed in Italy [19]. We assume that an average COVID-19 patient transmits infection to 3.5 other people and that some patients can transmit infection to a far greater number of people. There is a hypothesis about superspreaders, and superspreader events have been described [21]. Such a mechanism was assumed for a HIV epidemic of 1980, which was explained by actions of a Canadian flight attendant who transmitted infection to hundreds of people throughout North America, although the virus found its way into the United States several decades earlier. The first confirmed case of AIDS was detected in St. Louis in 1966. However, the epidemics did not develop to a rapid spreading phase during the Summer of Love of 1969, but progressed slowly for 11 years. Records of the 19th century describe several cases of typhoid spreaders who did not have typhoid symptoms, but caused mass infection.

Trigger effects and pulsed activities are difficult to include in classical models based on systems of common differential equations. However, pulsed scenarios are important to consider with due regard of cases where many people simultaneously contracted infection, their long-term consequences, and the effect of inertia in the epidemic dynamics. Football matches, such as an Atalanta vs. Valencia match, triggered an acceleration of an epidemic. One event created a group of infected spreader fans in Bergamo, thus provoking a Λ-shaped peak of the epidemic with a 12% mortality rate according to current estimates. The mortality rate, in addition to daily new cases, is a key characteristic to be used to compare different epidemics.

What is comparable can be compared, that is, countries that are similar in the availability of healthcare and report unbiased statistics. A dramatic difference in mortality rate between neighboring countries (e.g., 10% in Spain and 4% in Portugal) is highly difficult to attribute to differences in the quality of medical services. An explanation can be based on the immunological imprinting hypothesis. If a HxNx influenza virus is the first one to cause flu in childhood, than the person becomes forever resistant to all strains that are similar to the initial type, but is unprotected against HxNy viruses and becomes ill every year because the immune system fails to detect any difference and activates the memory B cells that remember the old antigen. As is known from the history of vaccinations, cowpox provided resistance to natural smallpox, and their similarity was thus advantageous to humans. A series of experiments in mice confirmed the effect of a blind spot in the immune system; i.e., there is a threshold in detecting variants of viral proteins [22]. Infection with a new influenza virus strain leads to production of antibodies against protein variants of the old virus strains. Antibodies that bind their target poorly interfere with recovery. Hemagglutinin/neuraminidase virus types were considered in the hypothesis. However, human coronaviruses have been studied far less comprehensively, coronavirus vaccines have not been tested, and uncertainty may act as a factor here to increase differentiation of disease scenarios. There is no guarantee that antibodies to a desired protein always occur at a maximal frequency in the total diversity of possible antibodies. Several β-coronaviruses have long been present among the agents of seasonal viral acute respiratory infections, but the distribution of coronavirus types was not analyzed in regions as thoroughly. The hypothesis suggests that people who did not have coronavirus infections earlier develop a milder COVID-19 disease compared with people who were infected earlier with other members of the family Orthocoronavirinae. Older people might have an infection with a strain that is similar in antigenic properties in the 1980s, and only a weak response is now produced by one of the components of their immune system. An irregular frequency distribution as a probabilistic component thus arises in the formation of the immune response upon antigen presentation. Transmembrane spike proteins vary substantially among RNA viruses. Protruding spike regions are a priory thought to be immunodominant, while conserved nucleoproteins and large proteins responsible for RNA packaging in virions are also immunogenic (targets accessible to T cells). Immune memory cells that remember nucleoprotein antigens can act as a helpful factor upon infection with a virus that has modified regions in crown spikes.

To check the hypothesis of a blind spot in immune memory, it is interesting to compare the antibody variants, including those to two parts of the spike protein and alternative N proteins, between patients with mild disease and those with a fatal outcome and to analyze their frequency distribution. If negative tests for antibodies are observed in a substantial portion of severe COVID-19 cases confirmed by PCR tests and positive tests for antibodies are observed in PCR test-negative cases, the discrepancy will provide indirect evidence for the hypothesis. An argument against the hypothesis of an “original antigenic sin” is that circulating virus types must have been regularly distributed throughout Europe in globalization conditions. In genetics, variants of histocompatibility gene regions that are responsible for immunity are analyzed in terms of frequency to explain the differences observed in COVID-19 mortality rate.

Several methods are used to model disease spreading, but they are not universally suitable for prediction in a pandemic and, even more importantly, estimating the future mortality rates. Modifications of the SIRS model for spread of disease are common and have been studied both in the form of systems of differential equations and in a discrete algorithmic form. The algorithms included both fully deterministic ones [22] and those allowing for probabilistic parameters of outcomes (infection, recovery, and death) [23]. There are cameral SIRS models, which include the isolation of subpopulations and various risk groups for the spreading of special infections. Infection with a sensitive variant and infection with a resistant variant were separated in a study of certain treatment aspects [24]. A patient has a probability to die or become immune in automatic SIRS models, but it remains unclear how a transition to the asymptomatic carrier state should be included in the transition set. Young people with mild symptoms remain virus spreaders (with viral RNA carriage confirmed) for 2 months according to reports from Italy; however, data on the time course of the immune response are necessary for modeling.

There are still no uniform guidelines on how restrictive measures should be introduced by a country’s authorities to stop the spreading of infection. Restrictive measures are efficient when introduced in all countries at the same time. The World Health Organization declared a pandemic on March 11, 2020, but after the fact. Predictions and mathematical models of coronavirus spreading scenario were available to the authorities, but did not play an appreciable role in decision making. The availability of numerous mathematical models [25] does not help when officers of the World Health Organization lack the will to make decisions and to introduce operating algorithms mandatory for all countries. Studies published as early as 2007 already warned that an epidemic was possible if a SARS virus circulating in chiropteran populations is transmitted to humans (see [26, 27]).

The issues of the optimal time to relax lockdown rules without provoking a repeated outbreak and when it becomes certain that a repeated increase in disease incidence is unavoidable in a given scenario are important to predict. Data from earlier epidemics show that repeated peaks occur 2–4 months after the first epidemic wave has waned, as was the case with Spanish flu and a plague epidemic that hit Marseille in 1720. Predicting a repeated wave of the epidemic is a pressing problem of modeling. It is important to distinguish a repeated outbreak as an event from an oscillation phase when the epidemic process follows an oscillatory scenario. The key issues that are of critical importance to know for modeling and prediction are how permanent immunity forms, how long it works, and whether a virus reactivation is possible when a latent infection site persists in the body. There are reports that the virus was detected again after its disappearance from the body in some cases. This is a concern, and it is unclear how rare such cases are.

One drawback of many SARS models is that qualitative modes of behavior are independent of the initial state observed at the start of a local epidemic process. According to one hypothesis, a difference in the initial number of infected people may explain the difference in epidemic scenario between different regions, acting together with other factors. The social contact density is the main one among these factors and is poorly formalizable in equations.

It is essential for modern epidemic modeling that the mortality rate is a functional rather than a parameter set a priory in a system of equations, even when it is set to be a random number with a known distribution through age groups. A model will have predictive potential when demographic and immunological parameters of subpopulations are used to estimate the probable mortality rate in computational experiments.

Detailed network models of virus transmission in a heterogeneous population can be built with current computational capacity, but the virus potential for latent spreading prevents the use of medical statistics. Conventional modeling methods are certainly of use for analyzing the local variants of how phases of the process and the portion of recovering patients change depending on the time course of daily new cases and the scenarios of the disease course in the body.

FORMALIZATION OF DELAYED REGULATION IN POPULATION PROCESSES

Phenomenological models are important when the actual mechanisms are unknown and infeasible to describe directly or when their complexity is impossible to formalize. The functions of immune subsystems were mathematically modeled in many studies [28]. Humoral and cell-mediated immune responses were long considered separately when constructing models. Marchuk et al. [29] analyzed immune system aging due to a higher antigenic load. However, the structure of innate immunity and difficulties in triggering adaptive immunity remained unclear until recently. Studies of the T-helper subtypes and receptor selection algorithms are still far from complete. The findings of Medzhitov et al. [30] mean that in terms of biological cybernetics that the system that regulates immune responses works in a trigger manner, through a cascade of signals, and belongs to the paradigm of discrete-continuous and hybrid types of dynamic processes [31].

In modeling the direct interaction between a pest and a forest or a virus and an organism, the problem is that the time scales of two fighting species are not comparable (in particular, a HIV virion lives for up to half a day). Assuming a classification of phenomenological models in biology, the interaction between the immune system and a virus does not belong to the predator–prey or parasite–host type, but is close to a model of war as a confrontation of two fighting armies. The efficiency of means of destroying the opposing side can be calculated in models without describing their internal mechanics. Scenarios of how the extreme interaction ends are also important to analyze in models.

Delay equations provide a common means to construct biological models with intricate dynamics [32] (a synonym of Volterra’s hereditary phenomena). In this approach, the regulation of an ongoing process corresponds to a past state observed in the system a certain period of time previously [33]. A delay reflects a set of actual time-related factors in biology [34] and is not restricted to the time of reaching maturity. As an example, several time-related parameters characterize various intracellular processes, such as protein expression mediated by mRNA [35].

Response delay factors are important in the case of formation of the immune response. The factors depend within a narrow time range on the frequency at which active cells encounter the virus. The body has both means to immediately respond to infection and means to produce a response that is specific to a particular target, but is delayed because dendritic cells have first to present the viral antigens (by analogy, as a composite portrait is presented by police).

Volterra was the first to consider the equations with N(t – τ) in mathematical biology. The Hutchinson, Nicholson, and Gopalsamy delay models are the best known in biology and have several modifications [36]. The Hutchinson equation was not derived by Hutchinson himself in 1948, in contrast to what may often be found in publications (old works have been digitalized and are available, and classic works are now possible to read rather than to refer to from habit). Hutchinson [37] briefly outlined the hypothesis that earlier states affect the reproduction efficiency; the idea was not the main idea of the work. The model was ascribed to Wright in [38], while Wright [39] provided a somewhat different form:

$$y'(t) = - \alpha y\left( {t - 1} \right)\left[ {1 + y\left( t \right)} \right],\,\,\,\,\alpha > 0.$$

May [40] was probably the first to write the equation in its conventional form:

$$\frac{{dN}}{{dt}} = rN(t)\left( {1 - \frac{{N(t - \tau )}}{K}} \right).$$

(1)

Equation (1) is actually a complication of the Verhulst logistic model with the delay τ introduced for the regulation that is determined by the carrying capacity K of the ecological niche. K is an important parameter and was theoretically grounded by the ecologist Hutchinson. K is more than merely a constant in a model, but is a principle of the intraspecific regulation in theoretical ecology [41] and is based on the a priori assumption that there exists the maximum allowable equilibrium population size $\forall t$. The mathematical model gave origin to common theories of r- and K‑strategies in species evolution, although a theoretical mathematical description of the European population of the mid-19th century was originally developed by Verhulst and was never intended for describing the extreme situations that arise upon new species invasion.

When the delay τ is low, the equation provides a phenomenological description for a mode actually observed for decaying population oscillations N(t) → K. At rτ ≈ π/2 ± ε, N(0) < K, trajectory (1) describes the character of a waveform epidemic, such as a repeated phase of growing new daily COVID cases in Texas (Fig. 4a), but not in New York or Italy. A distinct oscillatory scenario of the epidemics is observed on a plot constructed for Iran without averaging (Fig. 4b). The scenario can be called the “second wave,” but, mathematically, this wave is part of an integral decaying oscillatory process. Differences between the two situations were earlier described with the examples of a repeated pest outburst in a gum-tree forest as a separate stochastic event and continuous oscillations of a leaf-roller moth population in Canada as an integral long-term population phenomenon [42].

Fig. 4.

(a) A plot of new COVID cases in Texas with a smoothed curve obtained by data averaging over 7-day periods (data from The New York Times) and (b) an oscillating COVID epidemic in Iran.

The condition $\forall t$, N(0) < K, N(t) ≤ K is true for the equations included in the Verhulst, Gompertz, Richards, and other logistic models, but is inconsistent with rapid invasion scenarios. A scenario observed at $\exists {{t}_{\sigma }}$ < t_m, N(t_σ) > K, and N(t_m) < K corresponds to a classical σ-shaped curve, which E. Odum used in his Ecology to show a typical scenario where the size of an invader population grows with the overshoot N(t_σ) > K, but without an extreme outburst phase. The invader population cannot immediately stop increasing in size once niche K is reached (the population is unaware of this fact) and the overshoot N(t_σ) > K consequently arises until the size stops growing (N(t_σ = 0) and decreases without forced to an acceptable equilibrium with the biotic surrounding. The overshoot does not belong to extreme transition modes. The growth in size stops on its own (N(t_∞) = K), as a parameter regulated in the absence of external counteraction. The dynamics following Eq. (1) can describe the scenario only when an ecological niche is filled in a “peaceful” manner, without damaging the environment, via competitive displacement and the population size does not drop dramatically after the overshoot, that is, an increase over the a priori niche capacity. Cases are known were the invader population was too aggressive and died completely. The achievement of the point of a zero increase in population size does not necessarily mean a balance with the environment in many actual situations. As an example, K as the limit number of cells accessible to viral infection is determined by the efficiency of interferons in the body.

The solution of Eq. (1) is dissipative: N(t) ≤ (–1) + exp(rτ) regardless of N(0) [43]. The properties of Eq. (1) with the formation of an orbitally stable cycle and nonharmonic oscillations are well known [45] as the formation of autowaves in systems with a delay [44]. The trajectory of Eq. (1) passes to a cyclic mode for a population with the reproductive potential r after an Andronov–Hopf bifurcation when the bifurcation parameter product rτ > π/2 [46]. When biological problems are studied, Eq. (1) is analyzed with the constant prehistory $\forall t$ ∈ [–τ, 0], N(t) = N(0). A single-parameter form is used in mathematical studies: $\dot {N}$ = λ[1 – N(t – 1)] and the results are difficult to convert back and to express in terms of ecosystem processes. A generalization of the Hutchinson equation can be written as follows to allow for the age structure of the population:

$$\frac{{\partial N}}{{\partial t}} = r\left( {1 - \frac{1}{K}\int\limits_{ - {{\tau }_{1}}}^{ - {{\tau }_{2}}} {d{v}(s)N(t + s)} } \right)N.$$

(2)

The drawbacks of the Hutchinson equation and its generalizations are well known in applied ecology [47]. Relaxation cycle minimums (1) become prolonged and deep with the increasing amplitude and actually fall within the ε vicinity of zero. The mode is unrealistic in the context of interpreting the properties of an isolated population. There are grounds to state that the potential for essential biological interpretation of behavior is lost at high τ values so that the model fails to describe rapid abrupt changes.

An alternative blowfly equation does not explicitly includes the limit capacity of the niche:

$$\frac{{dN}}{{dt}} = rN(t - \tau )\exp ( - bN(t - \tau )) - \delta N(t).$$

(3)

The model in the form of a differential equation was proposed for describing the large-amplitude oscillations observed in a laboratory blowfly population in Nicholson’s experiments [48]. The well-known Ricker function is used in Eq. (3) to formalize intraspecific competition at an increasing density, but a delay is introduced at the reproductive term; i.e., it is an analog of the equation x_{n + 1} = $a{{x}_{n}}{{e}^{{ - b{{x}_{{n - 1}}}}}} - \delta {{x}_{n}}$. The Ricker function iteration that has long been used as a model of replenishment and harvesting of a fish population in fishery has a chaotic mode following the scenario of Feigenbaum’s period doubling. The iteration cycles Ψⁿ(x) = Ψ^{(n + p)}(x), p = 2^ϕ are unsuitable for describing the population fluctuations. A doubling of the period p creates pairs of points around the unstable point x*, while population cycles are serial permutations with a maximum at the end of a period: $x_{{p - i - 1}}^{*}$ < $x_{{p - i}}^{*},$ i = 0…p – 1. Analysis of Eq. (3) belongs to the field of chaotic dynamics [49] and is beyond the scope of this article. The Malthusian coefficient r already reflects the difference between the natural birth and natural death rates; then δ in Eq. (3) corresponds to additional mortality as an external factor.

Maximums arising in an oscillatory mode after a bifurcation of cycle creation in the Gopolsamy model [50] agree better with the views of regulatory oscillations:

$$\frac{{dN}}{{dt}} = rN(t)\left( {\frac{{K - N(t - \tau )}}{{K + jN(t - \tau )}}} \right).$$

(4)

Ecological models with a divergent argument were considered in many studies, and diverse modifications of delay models were reviewed [51, 52].

Model modifications are proposed below for describing the rare and specific biological processes that include variable scenarios. Unfortunately, epidemics do not occur in stable cycles. Oscillations of population processes are often sporadic and lack a constant period, thus representing intermittent fluctuations. The scenario of a spontaneous termination of a rapid invasion phase is important to describe in the case where reaching the equilibrium lim_{t → ∞}N(t) = K is an unreal scenario. Cyclic dynamics is not a main key property of epidemic models.

MODIFICATION OF A MODEL FOR CYCLIC ACTIVITY OF POPULATIONS

To obtain the oscillations that are realistic in terms of ecodynamics, a delay equation was develop on the basis of Eq. (3) with formalization of nonlinear threshold counteraction, which occurs in the mode of oscillations that arise after the Andronof–Hopf bifurcation and decay spontaneously. The equation is as follows:

$$\frac{{dN}}{{dt}} = rN(t - \tau )\exp ( - bN(t - \tau )) - q\frac{{N(t - 2\tau )}}{{A - N(t - \tau )}}.$$

(5)

Model (5) can describe the specific dynamics of the COVID epidemics in Michigan as a series of decaying peaks following a main abrupt outburst (Fig. 5). A similar development of the coronavirus epidemic was observed in Wyoming; however, the second peak was far higher than the first one. Such a scenario has long been known as sawtooth-like oscillations with different-sized decaying peaks in pest activity dynamics. Identical dynamics of phase changes can be absolutely different biological processes, which may be explained by a common mathematical form of their regulation functions [53].

Fig. 5.

(a) A series of decreasing peaks of the epidemics was observed in Michigan; (b) spontaneously decaying oscillations arise in model (5) after a bifurcation because the parameter r increases.

The model scenario is important in the assumption that the epidemic will not wane in full, but will lose its pandemic character to become as common as other sporadically active virus acute respiratory infections. Pests do not disappear after their outburst is over; however, their population size becomes small and their activity is not included in forest reports. In model (5), an abrupt change in N(t) without a bifurcation can evoke a new transition mode with oscillations.

Comparisons confirm that decay in an epidemic after a peak is a rather volatile, nonmonotonic process with a stochastic component and local outbursts. A smooth decrease is suggestive of willful falsification of statistics for reports. The variance of a series decreases when statistics are falsified. A scatter of more than 25% in daily new cases is commonly observed in real epidemics.

To further modify the model, it is reasonable to hypothesize that a relative regulatory factor, rather than a constant factor with a value fixed a priori, is present in the ecological niche. A transition to an outburst cannot start with single individuals in natural insect populations. It is important to note that data on forest damage, rather than the pest population size, are shown in outburst dynamics plots. The reproduction efficiency is low and depends on stochastic events in small groups. A sufficiently dense group of individuals is necessary for triggering an outburst, which may result from their efficient reproduction. A simple modification of Eq. (4) with nonlinearity increasing in a nonuniform manner in the numerator and denominator upon the regulation of reproduction and a relative character of the niche is:

$$\frac{{dN}}{{dt}} = rN(t)\left( {\frac{{C - {{N}^{2}}(t - \tau )}}{{K + \gamma {{N}^{3}}(t - \tau )}}} \right),\,\,\,\,C < K.$$

(6)

Thus, Eq. (6) makes it possible to obtain a sawtooth-like oscillation mode (Fig. 6) where a series of population size peaks start after the Andronov–Hopf bifurcation in the vicinity of a steady state that exists before the bifurcation and where the first triangular peak of the transition mode is greater in amplitude than the subsequent peaks. To trigger an oscillating outburst, a primary inducing factor is necessary in the case of Eq. (6). A certain effect exerted on the regulatory mechanisms makes the equilibrium that existed in the system lose its stability.

Fig. 6.

Λ-shaped peaks of an oscillating outburst in model (6) without changes in parameters.

To allow for counteraction of the environment to invasion of an aggressive species, the additional nonlinear term Ψ with a delay is introduced in the right part of Eq. (6):

$$\frac{{dN}}{{dt}} = rN(t)\left( {\frac{{C - {{N}^{2}}(t - \tau )}}{{K + \gamma {{N}^{3}}(t - \tau )}}} \right) - \Psi \left[ {{{N}^{k}}(t - h)} \right].$$

(7)

Equation (7) makes it possible to obtain a series of different-sized peaks and oscillations with a nonstationary period. A specific mode corresponds to the scenario where an oscillating outburst with intermittent Λ-shaped peaks is suddenly triggered in an insect population after reaching a substantial average population size. Figure 7 shows the dynamics of forest damage area in two regions of Alaska (LACL—Lake Clark National Park, KATM—Katmai National Park). The activity of the pest beetle Dendroctonus rufipennis is responsible for the damage and has been observed over a decade to occur as an outburst with a series of Λ-shaped peaks. Canada and Alaska have reliable public information on the state of forests over a century and the extent of damage to forests by various insects. Canadian data are therefore better to use in verifying the adequacy of a model when modeling extreme activity in biological processes. The dynamics shown in Fig. 7 confirms the hypothesis that population oscillations observed for insect pests of forests in nature are unstable transition modes and sporadic events in ecosystems, in contrast to the stable cycles of Volterra’s systems of equations.

Fig. 7.

The dynamics of damage to forests in two Alaskan regions after the start of an oscillating outburst of the spruce beetle Dendroctonus rufipennis (from a Forest Service report [54]) https://www.nps.gov/articles/insectsswan.htm.

Another regularly oscillating activity (Fig. 8) is observed in the mosquito Aedes aegypti, which is a carrier of dengue fever in South Asia. The threshold at which a new peak starts was found to increase with warming [55].

Fig. 8.

Activity peaks of dengue fever carriers in southern China (from [55]).

In model (6), it is important how close the size of the invader population is to a certain subthreshold level at which the pressure of the invader (or virus) population on the environment becomes high. The level differs between different biosystems (or organisms). Outbursts with series of rapid peaks (sawtooth-like oscillations), which were observed in certain insects, belong to types of extreme dynamics. Scenarios of this kind are beyond the theory of density-dependent regulation and the tendency to equilibrium in biosystems. The relative equilibrium niche capacity becomes the lower value for outburst activity. Regular series of outbursts usually tend to decay in the butterflies Malacosoma disstria and Lymantria dispar. The similarity of regulatory functions in these phenomena is important here.

THE INTERPRETATION OF THE AGGREGATE DELAY

The delay parameter is independent of the characteristic of a biological species and the period durations in its life cycle. In New Brunswick, the spruce budworm moth Choristoneura fumiferana generates extended activity cycles [56] with decade-long intervals between catastrophic outbursts. The intervals between peaks do not correlate with the duration of the development of an insect generation. The delay that is included in population models (1) and (2) is a phenomenological description of changes in population size dynamics [57]. This method extends the possibility of describing the changes in modes of trajectory behavior, but with changes in parameters. Known equations do not include the factor of changes in qualitative dynamics of scenarios as a function of the initial conditions [58]. The dynamics of infectious diseases do not necessarily follow one scenario because the immune system is individual to each person and there are unlikely variants and particular cases.

The delay τ, rather than the reproductive potential r, is a bifurcation parameter in Eq. (6) in my opinion. A delayed regulation is an adaptive aggregate characteristic in the model. The delay τ results from the interaction that forms between the species and particular biotic environment in the given region. The acting delay value in Eq. (5) includes the ontogenetic delay of developmental stages, the resource restoration rate, the time to response from parasites and pathogens, and migration. It is expedient to isolate the regulatory reproductive and adaptive components in the delay when analyzing the models. The latter is necessary for the population (or the organism) to develop adaptive responses and may change during invasion. In the case of a confrontation between the immune system and a virus in the body, the delay is an individual-dependent characteristic of all activation reactions and its average value may vary in local subpopulations (not necessarily as a result of genetic factors).

The Hutchinson and Nicholson equations fail to describe the end of an outburst or an epidemic, which cannot be a stable cycle and continue at t → ∞ in reality. Outbursts of insect populations and virus epidemics are transition modes, which are limited in time and have phases in their development and termination. Model (5) shows a possible termination path, but lacks variation. All extreme biological processes inevitably end, even if with complete degradation of the environment and death of the intensely reproducing species. Several scenarios are possible for termination and, in particular, the situation where minimal equilibrium is achieved with the formation of a latent infection site that does not exert an appreciable effect on the environment. As was noted previously in a study of a trigger iteration model of population outbursts in psyllids [1], describing the scenario of triggering the eruptive phase and modeling the end of an outbursts are two separate problems, and a spontaneous decrease in population size is more difficult to model. A virus imposes less strict requirements on its environment, and activation of virus suppression mechanisms will exert a triggering effect in the case of infection. We consider a bifurcation-free model of infection development scenarios.

MODELING OF DIFFERENCES IN SCENARIOS OF THE INFECTION PROCESS

A complex system of interacting cells of many types governs the set of body responses to suppress infection. The immune response is a cooperative response. A narrow role, or a few roles, as in the case of B cells, can be defined for each component [59]. Certain steps must be performed in the activation algorithms of each component of the immune response, and certain time is necessary for their activation. There are two large interconnected structural complexes, innate and adaptive, in the immune system. Different times to activation are the main difference between them in the context of modeling the infection scenarios. Antigens should be captured and presented to helper T cells (the distinguishing characteristics of the antigens should be detected), and respective antigen-specific killer cells should start propagating. The processes cannot occur in an instant. Antigen presentation is a process that results in migration of mature T cells into tissues to recognize and kill the infected cells that express the respective antigen on their membranes. It is problematic to reduce the time to a response from the adaptive cell complex.

Virus replication is a continuous process in contrast to a succession of insect generations, and a delay differential equation is therefore proper to use in order to analyze the scenarios. The succession of phases in the course of a disease is interesting to model based on the internal logics of counteracting invasion in scenarios that lack external induction and bifurcations due to fluctuations of internal parameters. Changes in the parameters are possible in the long term. It is irrational to describe transitions between phases in primary infection by fluctuations in characteristics and apparent bifurcation-related changes.

The main problem to analyze with the modified model is the issue of what scenarios are possible for spontaneous termination of an extreme process, that is, its transition to another quality after the phase of a rapid growth in virion number. The initial number of virions entering the body is thought to be important in certain situations [60]. The disease course may depend nonlinearly on the initial N(0) [61]. A dose-dependent immune response should be considered in many cases, as is important for efficient vaccination [62]. Dendritic cells have to interact with lymphocytes after contact with an antigen. When N(0) is low, the virus may evade the immune system.

Two (low and extremely high) doses of influenza virus were used to infect different groups of mice in experiments [63]. A 1000-fold difference is of interest in experiments, while in reality a tenfold difference is already large. The initial infectious dose was shown to regulate the early response and the release of immunomediator cytokines. At the higher initial dose of virus, symptoms of inflammation were more intense, interferons were induced more rapidly, and, unexpectedly, active infiltration with CD8-positive T cells was observed earlier in the lungs. It was noted that, surprisingly, although virus titers in the lungs were substantially higher in the mice infected at the high dose in first days after infection, relative virus amplification at a peak viral load was substantially (125 times) higher in the mice infected with the lower virus dose. The initial infectious dose affected only the response from cell-mediated immunity according to the experimental findings, while production of neutralizing antibodies remained independent in its dynamics. Surviving mice were infected again, and those that had survived the higher virus dose were found to acquire resistance due to persisting CD8 cells.

Unfortunately, integral parameters of the immune response were not used in the study, while a generalized plot would make it possible to analyze the problem mathematically and to build the functional dependence between the initial infectious dose and the response from the immune system. The response should be similarly suppressed later because IgG4-class antibodies decrease inflammation by inhibiting the activity of phagocytic cells [64]. Threshold effects are the most interesting in the hypothetical dependence between the dose and response rate. The initial infectious dose can determine the qualitative characteristics and rate of the immune response in virus diseases according to experimental data and should be considered as a factor that affects the disease scenario in model computational experiments.

Changes in the behavior of Eqs. (5) and (6) are observed when a bifurcation is induced by instantaneous changes in control parameters; the changes are not always possible to attribute to essential factors. It is necessary in each case to theoretically explain the mechanisms by which an increase in r (or τ) arises to provoke a fluctuating outburst. All characteristics of species are variable in nature; however, their changes are usually smooth.

A new model is proposed that assumes an important role for the subcritical pathogen level H, which is the subthreshold viral load: N(0) < H < K. A logarithmic, rather than quadratic or exponential, function is used as a regulatory function:

$$\frac{{dN}}{{dt}} = rN(t)\ln \left( {\frac{K}{{N(t - \tau )}}} \right)(H - N(t - h)).$$

(8)

In Eq. (8), H reflects a particular state of the body upon the interaction between the immune system and the given infection and K is a threshold set a priori, but is not the same as the carrying capacity of the niche in the logistic model. The limit saturation level is understood as a state that is achievable and may be slightly exceeded for a short period of time, but cannot be preserved in the environment. The characteristic plays a dual role, i.e., the lower its value is, the milder the symptoms are, but the closer the critical threshold is for which a transition over it is irreversible.

We consider two main scenarios of the new model (8). We assume that h can only slightly exceed the regulatory lag τ, i.e., h = τ + ξ; τ is a substantial value; and H = K/3 is constant. In the first scenario, the phase of a decelerating logistic growth in population size is identical to the Hutchinson model at rτ ≈ π/2 with the overshoot N(t_m) > K, N(t) → K, but equilibrium achieved after decaying oscillations is unstable. There is no final equilibrium in this scenario (Fig. 9).

Fig. 9.

The nondissipative trajectory by which the solution of model (8) passes from the state N(t) ≈ K into the mode N(t) → ∞ at τ = 58, K = 15 000, H = 5000, r = 0.00000335, N(0) = 190, and ξ = 2.

The trajectory t > t_m abruptly leaves the vicinity of K, so that lim_{t → ∞}N(t) = ∞ and the solution of model (8) at the given values of N(0) and τ is no longer dissipative. This scenario means destruction of the environment in the case where preexisting initial immunity is insufficient for controlling invasion. An acute viremia phase and the stage with t_K > 2τ that occurs after max  N(t) and consists in long-term pseudostable severe disease are followed by irreversible destruction of the environment. In a computational experiment, computation of the transition mode after a trajectory outbreak stops with an overrun error from the accuracy of the numerical method.

A different scenario is observed in model (8) where the initial population size upon infection is closer to H. After a similar growth phase with N(t) → K, there is a transition from short-term fluctuations in the vicinity of K to fluctuations that occur with a far lower population number in the vicinity of H and decay slowly (lim_{t → ∞}N(t) = H). The growth phase is more rapid, but the time of infection in the vicinity of K is shorter. In the computational experiment shown in Fig. 10, the trajectory reaches stable equilibrium after a transition mode with N(t) ≈ K at similar parameters of model (8). The condition N(0) > H was not considered because it is impossible in reality; certainly lethal virus doses are injected only to mice in experiments.

Fig. 10.

The equilibrium that the trajectory of model (8) reaches after N(t) ≈ K in the state N(t) → H; N(0) = 1099.

Alternative scenarios of asymptotic behavior are possible for the solution of Eq. (8) without changes in key parameters of the interaction. In the context of viral infection, r, τ, and H are constants. The initial population size N(0) is important in comparison with the subcritical capacity H of the environment. It is assumed that N(0) ≥ μ₀; i.e., a near-zero population size does not start the infection. The parameters H and K are interpreted as characteristics of the body, namely, properties of its innate immunity. K varies within the host population and its variation is responsible for the variation in severity of symptoms. H is more conserved, but fails to remain constant with age and is a natural variable environmental parameter that influences the behavior of the model. When the lag is reduced to 3/4 h in model (8), the resulting scenario has a less distinct acute phase of infection, that is, milder disease occurs (Fig. 11). Direct changes in parameters more substantially change the position of max N(t) in the transition mode.

Fig. 11.

The scenario of asymptomatic chronification of infection with max N(t) < K is observed at a lower regulatory lag in model (8); N(0) = 1200.

Symptoms seem to be insignificant and of a short duration in the transition dynamics because the stage with N(t) ≈ K is absent; however, the asymptotic state of the process is similar to that in the scenario with severe symptoms. Biologically, the scenario can be explained as follows. Innate immunity and antibodies quickly stopped virus propagation, but did not kill all infected cells, and long-term equilibrium was thus created. It can be assumed that the immune system became “blind” and failed to complete its function (without T-killer activity), but a more logical assumption is that the immune response slowed. Negative feedback is known for the immune reaction and is due to suppressor T cells, which inhibit activities of other lymphocytes. The third scenario is deceptive because patients may be unaware of the severity of their illness for a long time. A one-way increase in reproductive activity to r_m ⪢ r in model (8) logically triggers the fast mode with N(t) → ∞ as an irreversible fatal variant of development with a short period of staying in the vicinity of an unstable equilibrium.

The qualitative pattern of disease dynamics obtained in model (8) is most common for hepatitis C virus (HCV) infection. Figure 12 shows a transition from an acute to a chronic phase, which is observed in reality [65].

Fig. 12.

The dynamics of the HCV RNA concentration and generation of various immune response types as acute infection changes to a chronic scenario (according to [65]).

Plots constructed from the start of infection show how quickly alanine aminotransferase becomes detectable once the virus reaches the threshold vicinity of K. Alanine aminotransferase is an intracellular enzyme and is detected in the blood at substantial concentrations when liver cells are damaged; characteristic acute symptoms develop at the same time. A peak of alanine aminotransferase release is delayed relative to the growth of HCV RNA concentration in the blood. Interferon-stimulated mechanisms are activated faster [66], but are incapable of suppressing the virus and act as a damping regulator. There is a substantial lag in the appearance of efficient lymphocytes that are capable of decreasing the HCV virions and subsequently reduce the alanine aminotransferase level. Even a greater lag is observed for the appearance of specific neutralizing antibodies in this scenario. Protein expression from a large group of interferon-responsive genes involved in interferon-stimulated defense mechanisms (not all of the genes have been studied) suppresses virus propagation in cells [67], thus creating a threshold K level, that is, a level of temporary unstable balance where virus reproduction remains reliably controlled.

An increase in virus particles over the unstable level K is allowable for a period of time, as is seen in the scenario of chronification with decaying oscillations, but the period should not exceed τ.

Mathematically, intracellular antivirus proteins involved in interferon-stimulated defense mechanisms and the zinc-finger antiviral protein (ZAP) convert the exponential growth to a logarithmic curve. The regulators provide time for cell-dependent immunity to execute the response algorithm. Primary defense mechanisms always work from birth and would have saved the victim from infection, but certain viruses have evolved means to counteract interferon synthesis and to evade ZAP. Innate immune-system cells, which do not undergo antigen-dependent selection, fight many viruses poorly and prevent infection sequels poorly. A CD8-positive T-cell response produced in a timely manner is crucial for stopping infection; a flu lasts 1 week because the maximum T-cell concentration is observed 7 days after infection in experiments [68].

The counteracting response included in the above scenarios does not completely suppress the virus, but converts infection to a state where a persisting infection focus does not threatens the host’s life, although presenting a latent threat. The Ebola virus is hidden in Sertoli cells. The time interval between contraction of infection and the start of activity of specific immune system cells is an important factor in resistance. The initial dose of the infectious agent is another important parameter in model (8). The two delay parameters of the model are determined, first, by the properties of virus replication and virion assembly algorithms and, second, by the efficiency of a complex set of immune response mediators, that is, the interleukin system of mutual cell stimulation. A lag is an attribute of the replication–suppression system.

An equilibrium state does not preclude repeated peaks of infection. Many viral infections recur after prolonged immunosuppression. The immune response shown in Fig. 12 usually decreases with the decreasing HCV load, and only IgG remains in the plasma. To cause a disturbance in model (8), changes in the key parameter r may be unnecessary, and abrupt changes in H or τ suffice to cause a new transition mode with an infection growth phase.

VARIANTS OF MODIFICATION OF THE INFECTION SCENARIOS WITH THERAPY

Chronification in model (8) is an unenforced transition of the trajectory through the vicinity of an unstable state to stable equilibrium. In one important scenario the aggressive invader species dies in the environment by being destroyed by the species. An eradication model with N(t_∞) = 0 can be constructed using the idea of Bazykin’s equation, that is, by including the population size L that is clearly critical for the existence of the species, but raised to the power of 1/3:

$$\begin{gathered} \frac{{dN}}{{dt}} = r\ln \left( {\frac{K}{{N(t - \tau )}}} \right)N(t)\sqrt[3]{{N(t - \kappa ) - L}}, \\ 0 < \kappa < \tau . \\ \end{gathered} $$

(9)

Models with the critical lower limit population size L require that N(0) > L, and a complete elimination of the disease focus is obtained in this scenario (9), while variants of extremely fast invasion of an initially small amount of a pathogenic organism escape analysis [69]. Viruses do not need partners to reproduce and only have to evade the immune system. A single virus is capable of causing smallpox according to published estimates [70]. Fortunately, this is not true for common seasonal viral acute respiratory infections and COVID.

In the context of computational modeling, subsystems of the immune system are conveniently classified into three categories by time to response. A third category includes killer cells and neutralizing immunoglobulins with affinity for constant virus antigens. It is not necessary to know their physiological mechanisms for mathematical formalization, but it is important to know the time to their triggering and the corrective factors that adjust their effects. Quantitative parameters of the adaptive component of the immune response may be substantial, but the lag in their achievement may be critical for the body, leading to a fatal outcome. Insufficient specificity is observed in some cases; that is, the interaction with a given virus subtype is poor. To reflect active counteraction by means of special therapies and, for example, cytotoxic lymphocytes and CD8-positive cells, the nonlinear component Ψ[N^k(t – h)] can be expediently included in an infection model as follows:

$$\begin{gathered} \frac{{dN}}{{dt}} = rN(t)\ln \left( {\frac{K}{{N(t - \tau )}}} \right)(H - N(t - h)) \\ - \,\,q\frac{{N(t - 2h)}}{{{{H}_{1}} - N(t - 2\tau )}}. \\ \end{gathered} $$

(10)

Special cells are capable of eliminating intracellular parasites, which persist in sites that are inaccessible to antibodies. In addition, model (10) makes it possible to obtain the scenario where the asymptotic tendency N(t) → 0 is very slow. Lymphocyte activity starts to be inhibited by the body in reality [71]. It is interesting to consider the model scenario where the activity duration of immune killer cells, which undergo an intricate selection procedure and are active towards a particular target antigen, directly depends on the initial state at contraction. The component Ψ[N(t – h, N(0)] should be formalized for the purpose. Its formalization will be considered in a future publication, when data on coronavirus infection phases and the dynamics of immune responses of different types are available. There are all grounds to believe that a great diversity of disease trajectories will be reported and that several dynamic scenarios will be described for confrontation with the components of the immune system. Supporting data are accumulating for the observation that an imbalanced immune response with a dysregulation of the expression profile of inflammatory cytokines (not only for interleukin 6) [72] makes COVID sequels more severe.

It is expedient to model antiviral therapy differently, by using the functional r(t) to reflect a decrease in reproductive potential because the available drugs used to treat HCV and HIV infections suppress virus replication and inhibit virus polymerase, but without causing lysis. The method is important for predicting the dynamics after the end of treatment.

We note an important applied conclusion from the model scenarios: equilibrium is an asymptotically stable state in the chronification scenario, while trivial equilibrium is not. If therapy is started when the system is in a stable state, its efficacy is low and ultimate success is problematic. Therapy should start in a transition mode with an unstable balance and a cure may rapidly be achieved in this case without the disease becoming chronic. Mathematically, one possible solution is to convert the process to an acute phase again in order to trigger the immune response simultaneously with the start of therapy, but this solution will not be approved. An argument against the solution is that the state K of the critical unstable balance may be ephemeral when a repeated acute phase occurs in a body that is already exhausted by fighting the disease. In our case, the higher the intermediate value of N(t) is, the greater the symptoms and risks are.

CONCLUSIONS

Comparisons were performed for original models (5), (6), and (8) of specific trajectories of biological processes in biosystems. Different variants are possible for the start (and the end) of extreme or fluctuating dynamics of invasion according to the respective equations. Variants arise, first, when the reproductive activity parameter increases and bifurcations occur and, second, when abrupt changes occur in the invader population size without changes in other parameters. Third, transformations in properties of a transition mode depend on the initial infectious dose of the virus and its new environment in Eq. (8), which is important for analyzing the infection phases in the case where a lag is observed in the immune response. The infectious dose was varied fivefold in the computational scenarios. Bifurcation-dependent scenarios of changes in behavior are also possible in model (8). Transformations of limiting sets of trajectories are not as important for analyzing the rapid phases of invasion. Separate consideration was given to the models that allow for apparent nonlinear counteraction, which is described using an additional functional of the effect (variant (7)).

Alternative imitation scenarios were consecutively analyzed with new model (8), including development of the infection process to a fatal outcome, a critically severe course of acute infection with persistence of a chronic focus, and an asymptomatic scenario with max N(t) ≈ 2/3K. Two variants substantially differ in the transition mode and are similar in asymptotic behavior. The equation is relatively simple and easy to understand; there is no need to adjust the reproductive coefficient r of the virus, change other parameters, or induce bifurcations for each of the three scenarios. Computational experiments based on the formalization of a delayed interaction in the equation were performed by varying the initial conditions, which existed prior to infection, relative to the concentration that is subcritical for the virus and specific to the current state of the infected organism. The three scenarios do not mean that changes in the mode of behavior cannot be triggered directly by disturbing the coefficients r and τ in the model, but it is not as interesting to describe the transitions between disease stages by fluctuations in parameters. The structure allows a range of modifications for analyzing the scenarios of virus elimination with N(t) → 0 with an additional nonlinear active influence in model variant (10).

The parameter K included in model (8) to consider the infection process is not identical in meaning to the carrying capacity of an ecological niche or the amount of cells accessible to infection because a fatal outcome would otherwise occur far earlier than this level is achieved. The level K is the infection level that can occur in the system only for a limited period of time, but cannot be tolerated by the environment; this level is not constant in a population. The process consequently either passes to a terminal stage or decreases to a stable equilibrium level. The critical state K is not achieved when the immune system responds efficiently and quickly. It is clear that the mode N(t) → 0 will not occur in real processes. The variant where dissipative properties are lost by the model solution is considered to be a terminal scenario at N(t) = 5H. A simple and rapid fatal scenario, like N(t) = N(0)e^rt, is certainly possible in model (8) when the reproductive parameter r is high and not compensated for. The higher the values of K and r are, the more difficult it is for the body to overcome the acute infection phase. The time interval t_K, in which the trajectory occurs in the vicinity of the limit infection level, is an important characteristic calculated in model (8) for the transition process. The parameters t – τ and t – h, which were included in model (8), are aggregate characteristics.

In a previous iteration model [1], pseudostabilization occurs before an insect outburst at a threshold state. In the new model, a threshold state is stable for a long period of time, but equilibrium achieved at higher values of the population size K is repeller equilibrium. There is no extended period of a stably high population size in the case of a dangerous pathogenic virus, but a chronic focus may persist. Various scenarios of phase succession are possible for many infectious diseases even within a subpopulation, as is clearly seen with COVID-19. The properties of the model scenarios explain the observations made in Italy at the same average viral load. One mathematical explanation is that that limit critical saturation level differs, while the subthreshold asymptotic level is the same and is preserved for a long period of time. Average values were compared in [18], but data averaging fails to fully characterize the essence of a phenomenon in some cases, and the scatter of values within the groups to be compared is important.

Abundant data are now available for the dynamics of the immune response and the succession of phases in coronavirus infection, but the data are discrepant. Logic suggests that many diverse variants will be found to enhance the response. Immunosenescence is observed in elderly people with a lower interleukin-7 level and a lower sensitivity to interferon γ; i.e., T-cell proliferation stops and the immune response becomes insufficient. However, fatal outcomes were observed in healthy young people without concomitant disorders, many of them were healthcare providers. Deaths of healthcare providers are still the main argument that the infectious dose exerts an important effect.

There are reports that a hyperactive response of immune cells is associated with the most severe sequels of coronavirus infection for the body [74]. Cytokine release syndrome occurs as prolonged systemic inflammation, while a cytokine storm after a cascade release of signaling molecules is an irreversible extreme phase [75]. Normal cells also undergo apoptosis after fast hypercytokinemia.

We assume that the regulatory algorithm makes it possible to adjust the targeted response, be it rapid and massive or delayed. Slow response regulation, which the immune system utilizes in asymptomatic disease, is a rational and cautious method in this context. An affinity maturation mechanism is triggered and, algorithmically, consists in a spatial search for antibodies with higher antigen binding constants among the available amino acid sequences. In terms of biological cybernetics, when extreme features are lacking and the situation allows, an optimal strategy is taking time to produce agents with high affinity for antigens (that is, IgG3 antibodies and mature CD8-positive killer cells) rather than immediately releasing huge amounts of neutrophils and macrophages in the lungs. SARS-CoV-2 is most likely capable of entering these cells and killing leukocytes and thus triggers a destructive inflammation feedback loop in alveolar sacs of the lungs. The function of an intermediate component is necessary for the flexible regulation scenario. The role may be played efficiently by NK cells, which account for less than 10% of all lymphocytes. Initial amounts of means to control the response (IgG4 and suppressor T cells, which account for ≈5% of the total CD4 cell population) are important, but the respective tests are not included in conventional blood counts.

When the humoral immune response and antibody titers correlate poorly with the increase in viral load, antibody-dependent enhancement of infection it is likely; i.e., imperfect antibodies form unstable complexes with the virus and thus help the virus to enter leukocytes, be released, and start replication in monocytes and macrophages. The effect was responsible for the failure of a dengue fever vaccination campaign in Philippines. Disease symptoms in unvaccinated people without antibodies were milder than in vaccinated people with low-affinity IgG antibodies, which bind the virus poorly. Successful development of a COVID vaccine was not reported in China by the early summer of 2020, and it is most logical to assume that developers observed signs of antibody-dependent enhancement of infection (which increases the risk even when observed in 5% of cases). Modern achievements in molecular 3D modeling will help to solve the problem by selecting the optimal protein region and constructing a gene for its expression in a vaccine strain.

There are examples where severe sequels were avoided by controlling the immune response in a timely manner [77]. The example of an efficient response control that attracted the greatest attention pertains to simian immunodeficiency virus, which replicates without causing immunodeficiency in certain species (e.g., the mangabey Cercocebus atys), but causes symptoms similar to human AIDS in other primates. Immune activation was found to be limited in resistant monkeys, and their CD4-positive immune cells did not start proliferating. Their dendritic cells showed a low cytokine release after virus administration in experiments [78]. Mathematically, this is the following scenario in the notation of the models: $\forall N(0)$ < K: lim_{t → ∞}N(t) = H, $\forall t$ > 0, N(t) ≤ H.

Young monkeys developed an optimization solution for the problem of simian immunodeficiency viral infection without fighting the virus. With optimal control of immune responses, we will have a perfect means to overcome new threats without an acute phase. The SARS-CoV-2 genome [79] is large compared with the genomes of related viruses and codes for several new short proteins with unknown functions [80]. Artificial expression has been achieved for 26 out of its 29 proteins. Five viral proteins serve to overcome the immune system at the intracellular level [81]. Physically, RNA replication and packaging will take a longer period of time and inevitably proceed with stochastic errors to produce half-functional virion variants in the body, thus increasing antigen differentiation of infection.

It is interesting to experimentally study to the issue of to what extent the initial infectious dose may increase the risk of a fatal autoimmune scenario. As already mentioned, cytokine induction in mice infected at a higher dose was far greater than in mice infected at a lower dose of influenza virus [63]. If a strong relationship is assumed for the initial virus dose and the (maximal) concentration of main cytokine types, a step, rather than smooth, function will reflect the relationship. However, a smooth response is observed normally [82].

It is noteworthy that an initial lag in cytokine expression was observed prior to hypercytokinemia in experiments with tularemia [83]. A threshold-dependent triggering is therefore the most likely for the process. A storm scenario is mathematically possible when a threshold is reached suddenly and natural immunosuppressive mechanisms fail to control the fast chain activation of all other cell types. A trigger-based description of changes in the development of scenarios is proper to use when transitions are extremely fast. In the case of COVID, the initial dose is not only important, but it may act as a critical risk factor for a fatal scenario mediated by rapid hyperactivation of inflammatory reactions by the host immune system. A threshold value is often closer than expected.

The model with achievable unstable equilibrium K is proper for applied problems of ecodynamics, such as a biological control of dangerous invader species [84], fish resource management [85], and biological remediation of water bodies after pollution with oil [86]. The equation provides broad opportunities of further modification in order to evaluate the strategies of antiviral therapy and the infection scenarios with due regard to vaccination, which reduces the time to efficient response from cell-mediated immunity. The majority of modern vaccines are evaluated by titers of antigen-specific antibodies in the serum or on mucous membranes after vaccination. However, the response rate of immune memory cells was found to be more important [87]. A bottleneck was noted in the function of the most efficient component of the immune system; i.e., the rate of T-cell migration to lymph nodes may be limited by the throughput capacity of venules with a high endothelium [88].