key: cord-1048757-cepl51m8 authors: Perevaryukha, A. Yu. title: Simulation of Scenarios of a Deep Population Crisis in a Rapidly Growing Population date: 2022-02-10 journal: Biophysics (Oxf) DOI: 10.1134/s0006350921060130 sha: 37651804c4de6211eb509fa0af6cad7f62575d49 doc_id: 1048757 cord_uid: cepl51m8 Abstract—This article focuses on the modeling of crisis and threshold development of the population process during the formation of a new population in a competitive environment. As a population spreads, a deep population crisis may arise as a result an abrupt triggering of biotic countermeasures before resources for a further increase in population size are exhausted. A bottleneck occurred in the history of many populations, including humans at the time of the Neolithic crash in Europe. Invaders with high reproductive potential often exert deleterious effects on biosystems. The emergence of efficient competition can not only cause classical cyclical fluctuations, but also lead to a complete extinction of the population after a series of high peaks in its abundance. Two alternative scenarios provide classical examples of induced population crises. One was observed in Gause’s experiments where an introduction of a predatory ciliate drove another ciliate species to extinction. The other scenario was observed in a series of experiments where bacteriophages were introduced into colonies of actively dividing bacteria that had a dynamically adapting antiviral mechanism. In this work, modifications to the model were proposed to describe the actual scenarios of crisis effects in population dynamics. Equations with deviating arguments in the time variable allowed a threshold effect of conditions on reproduction of the invasive species and an aggregated nature of the lagging regulation with two time factors. The computational scenarios described both completion of the process after a repeated outbreak and successful elimination of the population crisis via rapid adaptation. Deep crisis phenomena are characteristic of local population dynamics when organisms interact with viruses that are new to them. This article continues works to improve the mathematical biological methods to analyze volatile and rapid transition processes. Recent studies have confirmed the majority of ideas about modeling the scenarios of virus infection spreading [1] . Many unobvious factors whose significance was hypothesized in early pandemics previously [1] have been established in detailed studies of the new disease or received experimental support by April 2021. Certain other questions that arose at the beginning of the pandemics remain unanswered. As an example, there is still no reliable explanation for why people of the Far East are, on average, more resistant to the virus than people of southern Europe, where the mortality rate is extremely high. The infection potential is still unclear for patients who have recovered from COVID, but remain PCR test-positive for a long period of time. Russian specialists have published several interesting studies [2] [3] [4] to contribute to investigating the important dynamic properties and prospects of the current epidemiological situation. The objective of this work was to mathematically describe the scenarios of a population crisis, which may spontaneously arise after a rapid growth phase even before the available ecological niche is filled by the species. The variation of the phenomenon is important to study for solving applied problems, such as designing methods to biologically control dangerous invaders. The principles of divergence in the dynamics of the COVID-19 epidemiological processes and its inevitable variability were noted in April 2020 on the basis of analogies. The role of superspreaders in triggering the COVID epidemic process is a feature of the current pandemics and has repeatedly been confirmed in various countries. In the context of modeling, the prob- lem is that many patients do not transmit infection to anybody, while some patients transmit it to many people. The superspreader effect has been known for many years and can be described by the example of the HIV epidemics that started in 1980. HIV found its way into the United States as early as the 1950s, but a single hyperactive spreader was enough to greatly accelerate the dynamics of the epidemic process in 1979 and 1980. Important findings are that superspreading is characteristic of coronavirus to a greater extent than of seasonal influenza and that the peak of the viral load often does not coincide with the peak of symptoms [5] . As a result, the number of COVID superspreaders is more important than the population-average basic reproduction number R 0 in a local case at an early stage. The dynamics of the early stage greatly varied among regional epidemics in Europe in early 2020 mostly because of the presence of active superspreader groups, which is hardly possible to foretell when predicting the scenarios. The second wave observed in the COVID epidemic cycle in the United Kingdom in winter 2021 (and presumably caused by the UK strain) was far greater than the first one in the number of cases, while the mortality rate was appreciably lower than in spring 2020 (Fig. 1) . The phenomenon is important for predicting the epidemic character. The decaying wave dynamics of the epidemics is qualitatively similar to pulsed sawtooth-like outbreaks of insect forest pests, whose dynamics has previously been described with an complexly customizable equation with four inclusions of various time lag values [6] : (1) Equation (1) generates decaying peak series wherein a second peak in a series can exceed the first one depending on N(0). As noted previously, the oscillatory character of the COVID-19 epidemics is necessary to distinguish from the scenario where a Λlike outbreak repeats and the distinction was confirmed in fact in Iran, New York, the United Kingdom, and the Republic of South Africa. The variant with a prolonged epidemic plateau in place of two waves was far rarer in local regions. A single rapid Λ-shaped epidemic outbreak was detected in Sri Lanka in 2020, being similar to the outbreak of H1N1 influenza in Philadelphia in 1918. The evolution of COVID proved to be more complex than that of Spanish flu. A classical pattern of a two-peak oscillating epidemic was observed in the Republic of South Africa in 2020 (Fig. 2) . The emergence of a new strain did not change the dynamic mode of oscillations and it can be predicted that a third peak of disease incidence will not exceed the first one. All large-scale pandemics end at some time, even with variable viruses that cause acute respiratory diseases. The immune response will not lose its efficiency due to CD4 and CD8 T cells in the period between peaks, resistance will accumulate 1 1 ( 2 ) ( )exp( ( )) , ( ) . A N t in the population, and the amplitude of oscillations in disease incidence will stabilize in the Republic of South Africa. It has been confirmed that T-cell memory and the effector function of CD8 cells efficiently work against many coronavirus epitopes [7] . Moreover, cross-reactive T cells have been observed in individuals without confirmed SARS-CoV-2 infection [7] . Emerging variants of the spike protein of the virus will certainly not all cause outbreaks of severe COVID cases. T cells of patients recovered without immunosuppression will recognize the epitopes of new strains, such as South African B.1.351, with a probability higher than 0.5. It is important to identify the situations where qualitative changes arise in epidemic dynamics, and hospitalization statistics are better for this analysis. The scenario that has a period of an extreme growth in daily new cases is far more dangerous and was observed in Brazil in March 2021. We note that a decay phase with a decrease in oscillation amplitude, which had peaks in August and January (as in the Republic of South Africa), abruptly and quickly changed to another mode in March, displaying a Λ-like outbreak with an exponential growth. The new outbreak apparently exceeds the initial disease spreading rate. Figure 3 shows the wave dynamics of hospitalizations in the southern Brazilian state of Rio Grande do Sul. Two weeks later, a similar epidemic outbreak started in Chilean cities after a quiet period in mid-March 2021. Brazilian strain B.1.1248 will probably be the most problematic in the so-called era of war between vaccines and strains [8] . The situation observed in South America suggests a trigger event that could abruptly change the character of the virus epidemic process 1 year after its start; this phenome-non must concern researchers. It is interesting to find the portion of reinfections among the cases observed in March 2021. The scenario of changes in epidemic dynamics (Fig. 3) is a bifurcation scenario and additionally displays signs of an abrupt loss of the established mode. This qualitative transformation of the pattern observed in southern Brazil suggests an abrupt increase in the basic reproduction parameter r. The increase was due to a mutation and occurred earlier than in early March, because the process is inertial. There is a lag in the response of hospitalization statistics. A model of the scenario that includes a similar disruption of the established cycle after the Andronov-Hopf bifurcation and the trajectory appearing in place of a new cycle at has been obtained with an extended modification of the delayed logistic equation [9] : (2) Cycle disruption in a nondissipative solution of Eq. (2) describes the moment when the process undergoes transformation because of an abrupt increase in r, but not the variants of its end. The scenario with a bifurcation is actually unpredictable and reflects an independent stochastic event. The H1N1 swine flu epidemic of 2009 occurred in these regions with two peaks, in July and November. Thus, an interval shorter than 2 months between a past peak and a new increase in disease incidence is a sign that the virus has acquired (accumulated several) dangerous changes. A typical dynamics with 6-month intervals between peaks is observed in Iceland: April, October, and April. Seasonal outbreaks of acute respiratory infections are not observed in Brazil or India, but are an important factor in Northern Europe and certainly Russia. Pooled case data are poorly suitable for analysis when originating from countries with large populations (the United States, Indonesia, India, and Brazil) and an uneven population density. Plots are actually combinations of dissimilar groups of local scenarios. An epidemic dynamics with three increasing waves is observed in the United States according to pooled data, but qualitatively different patterns are detectable in New York and Michigan. The methods to report statistics have been improved in Mexico, but a plot of daily new cases in Mexico qualitatively differs in dynamics from the respective plots of New York, the Republic of South Africa, and Brazil (Fig. 4) . The epidemic process developed in Mexico more slowly than in the adjacent countries, while the population density is extremely high in its capital. A comparison of various scenarios confirmed that lockdown measures fail to stop a rapid increase in disease incidence when this phase has already started and a lockdown, which is disastrous for the economy, is implemented only because there is risk to exhaust hospital bed capacities. This is an example of attempts at controlling a process on the basis of data that have a lag in reflecting its dynamics. Apart from lockdown measures, social norms, and communication traditions, particular immune response variants affect the dynamics of regional epidemics and have specific features in some local groups. Sets of cross-reactive B cells formed as a result of previous contacts with another endemic β-coronavirus were found in fact in some patients with severe symptoms [10] . Memory B cells (MBCs) account for a small, but stable, part of lymphocyte receptor repertoires and do not facilitate the neutralizing antibody response according to [10] . MBC activation failed to Combinations of four previously known coronaviruses circulated in various regions, differing in composition. This circumstance determined the apparent local characteristics in the percentage of severe COVID cases and the percentage of deaths in spring 2020 as a result of the blind spot phenomenon or the "original antigenic sin" effect [11] observed in the immune response. The properties of the effect are such that the effect will lose its significance in repeated epidemic waves of 2021, but is capable of reoccurring in the future because of mutations that accumulate in viral RNA. Antibodydependent enhancement (ADE) occurs, wherein the virus-antibody complex is not tight enough and the virus escapes endosomal degradation in leukocytes, but infects them instead [12] . ADE could be local and episodic in 2020, but is capable of becoming a common cyclic regional effect in patients who have recovered from COVID. MBCs that responded previously to a no longer significant epitope become a hindrance after a season of acute respiratory infections, producing suboptimal antibodies when encountering the slightly changed antigen. As an example, mass vaccination against the dengue virus failed because there are four different serotypes of the virus; after vaccination, dengue was more severe in some of the vaccinated people than in the control population on average [13] . ADE arises stochastically in patients infected with the dengue virus; however, a 10% failure rate after vaccination is an important problem. An accumulation of genetic differences and generation of new SARS-CoV-x strains were apparent as early as a pandemic was declared by the WHO. Coronaviruses possess a mechanism that provides the fidelity of their RNA replication, and such a virus may remain stable for many years, circulating in a local population. However, the law of large numbers plays against humans when the virus spreads through all continents. Deletions have come to contribute to virus variation; i.e., some regions are lost from the large virus genome. Unsurprisingly, the spike protein is affected by mutations, although there are conserved regions in this complex transmembrane protein as well [14] . I decided to get vaccinated with the understanding that there is a nonzero probability that another strain will circulate in January 2022, and IgG antibodies produced after vaccination may have insufficient affinity for a new variant of the receptor-binding domain (RBD) of the mutant spike protein. However, the larger the portion of virus-resistant people is in the global population and the quicker its accumulation, the lower the likelihood is of dangerous events, that is, new significant mutations acquired by the virus. The majority of events are not fully determined in the immune response of the body and the likelihood of success is still greater. Herd immunity resulting from mass vaccination can hardly vanish at once with the advent of new strains, but should decrease linearly as a result of similar stochastic factors that accompany the antigen presentation process. Immunodominant peptides of B and T cells will vary and have an complex distribution pattern in a large group of vaccinated individuals. Immunity will not lose its efficiency totally and at once. As was shown in many studies, it is important to investigate the cell-dependent response of CD8-positive T killers to a range of antigens that are regions of conserved structural proteins and nucleoproteins of the virus, but there is an inconstant time lag in the development of this response. A special role of natural killer (NK) cells in the immune response followed from the hypothesis of a lag in a T-cell activation scenario and was high in fact. A decrease in the effector functions of NK cells was observed in patients with severe inflammatory complications [15] . A comparison of new findings from several studies makes it possible to understand the reason that this cell group is important. Each component of immunity performs a range of specific functions, but their functions may overlap. NK cells provide antibody-dependent cytotoxicity, while alveolar macrophages are responsible for complications. IgG antibodies that target the S protein of SARS-CoV-2 and lack fucose residues in the Fc region (afucosylated immunoglobulins) have higher affinity for Fc receptors and consequently induce macrophages to release excess proinflammatory cytokines, which are signaling molecules [16] . Their release increases systemic inflammation via positive feedback. The greatest problem is to support the model hypothesis of a dose dependence in the scenario of an infection process, where the initial infectious dose N(0) strongly affects the development rate and quality of the cell-dependent immune response. The dose, but not a relative estimate of its virulence, can be used as a quantitative parameter in a computational model of a scenario. Direct experiments on humans with higher doses are impossible for obvious reasons. Certain data make it possible to believe that there is a threshold initial infectious dose that is capable of causing severe COVID by activating the cascade of inflammatory cytokines, including IL-1β, which trigger both local and systemic inflammations. Experiments with re-infection of rhesus monkeys showed that antibodies can protect one from SARS-CoV-2 infection in a dose-dependent manner, but only when the IgG level developed previously exceeds a threshold titer necessary for protection [17] . Experimental studies are necessary in order to evaluate the quality of the total spectrum of immune responses to various coronavirus doses, as in experiments described for various doses of a H3N2 influenza virus [18] . The developers of the VLA2001 classical vaccine, which contains the inactivated coronavirus, confirmed a pronounced dose-dependent immune response to their preparation. In the context of the complicated manner in which COVID depends on the initial virus dose, little attention is paid to using vaccines to improve the immune complex of mucous membranes, which pro-vide the first protective barrier and are important in individuals with compromised T-cell immunity and, consequently, increased COVID duration. Many separate findings, which are disaggregated and seem unrelated to each other, gradually form an integral mosaic of immune response variants, which are interesting to classify, in particular, using mathematical theories. The problem of the ecological dynamics of transition processes is as interesting mathematically, but has received insufficient attention. Spontaneous invasions have become regular or even inevitable events in the ecology of biosystems. Extreme forms of processes occur when imbalance arises in an established biological diversity and are of applied significance and interesting to describe mathematically. Simulation is described below for several alternative invasion scenarios that involve active invaders and are accompanied by an abrupt increase in competitive counteraction. Invasion scenarios differ from those of fights between the immune system and viruses in having a greater potential for mutual adaptation of the species, where the invader may win the competition for an ecological niche and become a dominating species. A virus lacks a stable "niche" and has only a minimal effect that does not destroy its environment. The body is incapable of existing for long in a state with a high viral load, as was reflected in a model scenario of the infection process exemplified by hepatitis C [1] . It is clear that infection of all available liver cells means a fatal outcome of the disease. Let an invasive species with a high reproduction parameter (r) intensely spread through a new region. At the first step, the dynamics of the process can be close to exponential with the growth equation N(t) = N(0)e rt -qN(t), q < r < 1. At a certain time point t = T, the species encounters active counteraction by autochtonous (or introduced for the purpose) biotic factors. There are only rare examples of situations where the invader species reproduces intensely and completely destroys its environment, as was the case with reindeer introduced to the Bering Island or humans on Easter Island. A trigger factor acts as a mechanism that regulates the species abundance with a time lag. Several variants are known for the development of such events, as well as interesting real examples. Human populations outside Africa can be considered invaders; they had to fight for their "niche." As was established in archeology, Homo sapiens populations that entered Europe and intensely colonized its northern regions experienced what is known as the Neolithic crash at approximately 6000 BC [19] . The long-term crisis lasted approximately 500 years ( Fig. 5 ) [19] . The causes and consequences of the long-term demographic decline are a matter of discussion in historical science. Equilibrium was established and persisted for a long time after the Neolithic crash, changing to a new trend of an increase in population size at approximately 4000 years BC. This is an example of a prolonged crisis, by historical standards, but not a very deep crisis. Such changes in abundance are known as boom and bust in ecodynamics. A mathematical scenario with a single short-term outbreak [20] is only its particular case. Changes with peaks and crises are rather typical and important because a bottleneck effect accompanies them. Genetic changes are associated with the effect because a small group of individuals gives origin to the population that develops after passing through a bottleneck. The number of allelic gene variants can be reduced to a dramatic extent, and such critical phenomena can be considered as evolutionary factors and causes of intraspecific differentiation. Various species were found to have a period with a small population size in their history based on genetic mtDNA studies. Dramatic fluctuations in population size are normal and exert no effect on fitness in other invader species. As an example, the warty comb jelly Mnemiopsis leidyi is a dangerous invader and uses cannibalism [21] to maintain its population when food resources are exhausted, thus being able to displace competitors while passing through population size minimums. The minimal population size L is characteristic of many, although not all, species. Neolithic human ancestors exemplify a favorable scenario with a prolonged minimum of N(t) → L, L > δ N and a low reproduction efficiency. Many species failed to overcome a sudden crisis, including invader Examples of population crashes were reviewed in [22] . Many of the respective populations had a high reproductive potential, as was the case with the ragweed leaf beetle Zygogramma suturalis introduced from America to the Stavropol region [23] . Experiments provided clear evidence for various scenarios of crisis dynamics. Gause's classical experiments with two ciliate species (1934) can be interpreted both as a situation where delayed counteraction is abruptly switched on and as a laboratory model that illustrates the behavior of an aggressive invader in a new environment. A tube with a Parameciuim caudatum colony that grew logistically (N(t) → K) was supplemented (t = t D ) with the aggressive predator Didinium nasutum, which is a natural enemy of the former. Elimination of the prey population was the final outcome of the interactions of the two ciliate species in experiments without refugia, and the predator population had a second population size maximum and then died in the absence of food (Fig. 6 ) [24] . The second maximum of Didinium population size was larger than the first one and coincided with a minimum prey population size. Several publications described the experiments as an experimental disproof of Volterra's model and the theory of a stable cyclic development of species in a two-level trophic interaction, but Gause's interpretation was somewhat different. Gause [24] noted incomplete determination for the laboratory experiments. That is, certain equations can be used to predict the course of the struggle for survival when the Paramecium and Didinium population sizes are large enough. However, various stochastic factors become highly significant when the population sizes are low at critical time points when one cycle gives way to another. As a result, it is impossible to determine the course of development in each particular case, and stochastic changes come into play again. Gause thus saw a phenomenon that had not been studied at his time and is known now as a stochastic blurring of the separatrix in mathematics. A certain disturbance in the initial conditions at N(0) ± ε [25] does not necessarily leads to a different asymptotic position of the trajectory when alternative attractors exist. An area with a set of scattered points, rather than the threshold L, forms the boundaries for domains of attraction. This property is observed in model (3), where H is the subthreshold virus load and N(0) < H < K for the given infection: The effect of separatrix blurring does not contradict experimental data and observations because an increase in the initial dose by n + 1 virion does not necessarily lead to a fatal outcome. Mice have a chance to survive even when infected with a huge influenza virus dose, as is known from experiments [8] . Equation (3) provides an illustrative, but not the only possible, example of stochastic properties that exist in a deterministic model. However, computational experiments with Eq. (3) cannot describe the spontaneous disruption of fluctuations in an emerging population cycle. Decaying oscillations and lim t→tK N(t) = K are observed in the first phase of a bifurcation-free scenario for t ∈ [0, t K ], but equilibrium K becomes unstable. At the time t → t K , N'(t) → 0 and N(t) → K, and the sign of the right part of Eq. (3) depends on the value E h = (H -N(th)), which acts as a disturbing factor depending on the trajectory state at the time point t K -h and can be positive or negative as deter- Mutual adaptation of the species is not implied in the scenario with introduction of a predator ciliate. The predator invader is counteracted mainly by its own food demands and uncompensated reproductive activity. The stochastic factors that Gause noted are associated with the sensitivity of his laboratory system to disturbances in its initial state at the time points t = 0 and t = t D . Experiments where the free prey dynamics is cyclic, rather than logistic, are of particular interest mathematically. Such experiments can be performed by adding the ciliate parasite Holospora undulata to the ciliate species. The crisis scenario under study is often associated with the start of an epizootic in the active invasion phase. There is another laboratory system in which species quite rapidly adapt to each other. A bacteriophage (virus) was added to an Escherichia coli colony growing according to the logistic law N(0), N(t) → K with a small overshoot N(t) > K and an unusual dynamics was observed [26] . The bacteriophage efficiently suppressed the growth of bacteria, but then the viability of the remaining bacteria increased abruptly. The system with the virus passed through a minimum and stabilized (Fig. 7) . The artificial system with E. coli-bacteriophage counteraction provides a laboratory model for real invasion situations. As an example, the gypsy moth Lymantria dispar was a harmful invader in North American forests and the pathogenic ∀ fungus Entomophaga maimaiga was introduced and greatly reduced the moth activity, but it took several years for the fungus to adapt to a new moth [27] . Unfavorable data were obtained in studies of crises because of E. coli adaptation. There was an intention to use viruses as antimicrobial agents in medicine in the 1930s [28] . It took several decades for researchers to understand the actual cause of the processes that occur in a tube with E. coli and the virus. The 2020 Nobel Prize was awarded for studies of the CRISPR/Cas9 mechanism, but in the field of genome editing rather than bacteriology. Regularly clustered interspaced short palindromic repeats (CRISPRs) found in the genomes of many bacteria provide virus signatures to the intracellular antivirus system that acts as molecular scissors and utilizes Cas endonuclease to recognize and cleave alien DNA. As is known to every user of a personal computer, signatures for antivirus tools should be updated on a regular basis. Horizontal gene transfer [29] allows bacteria to rapidly overcome an acute crisis and to establish a new equilibrium. Genome editing is now thought to offer unlimited opportunities from developing resistant potato cultivars to solving the problems of oncology [30] , but E. coli CRISPRs ruined the prospects of using bacteriophages in medicine. Thus, the two dynamically similar variants of a bottleneck crisis, which were considered above, are caused by different biological factors, but follow the common principle of delayed regulation and adaptation. In Gause's experiments, evolutionary adaptation was not implied by the scenarios, which implied only the generation of spatial heterogeneity. However, evolutionary adaptation is an important factor in the case of many introductions of insect and fish species [31] . The above crisis scenario underlies the development of certain adverse processes, and mathematical models are therefore interesting to obtain and compare for the two dynamics variants. The goal is to develop models that are capable of phenomenologically describe the extreme situations that occur after the first phase of active invasion. Invasions made the behavior of biosystems less predictable [32] . A harmless benthic invader can actually trigger a sequence of changes in the organic matter flow of a water body [33] . The history was previously described in detail for the idea of simulating with regard to both the time lag in biological systems and hereditary properties [34] , meaning that the regulatory functional of a process depends on several previous states of the system [35] . Four methods are used to obtain a computational description of the delayed aftereffect. The best grounded and most common one is based on differential equations with a divergent argument. Methods that are employed less often include Volterra's method PEREVARYUKHA with integro-differential systems and discrete functional iterations [36] with the recursive computational element x n + 1 = Ψ(x n , x n -1 ) -Q(x n -i ). A continuouseven representation of model time was proposed for hybrid models with regard for the effect that the stepwise character of ontogeny and delayed development at a lower individual growth rate exert on the viability of individuals [37] . Current achievements in developing modifications of continuous hereditary models in mathematical biology were reviewed in [38] . Certain features of using delayed equations are important to understand when providing essential biological interpretations to computational results, often in the form of relaxation fluctuations [39] . To choose the simulation method, it is necessary to understand beforehand that complications introduced in a model reduce the range of its parameter values in which the trajectory behavior is interesting and is biologically grounded. Population models with a large number of bifurcation parameters often demonstrate the nonlinear phenomena that are beyond the possibilities of essential ecological interpretation. Such effects with transformations of attractors, the appearance of "periodicity windows," and a doubling of the cycle period are associated with a transition to chaotic trajectory behavior, which is difficult to predict. Delayed equations were developed in mathematical biology to describe population size fluctuations, that is, cyclic regimens of various shapes [40] . Population cycles are known to form in isolated populations with high reproductive potential even in constant laboratory conditions [41] . In laboratory experiments, the fluctuation amplitude may be affected by the food replenishment rate [42] , as was the case with wellknown Nicholson's experiments with the fly Lucilia cuprina [43] . Three delayed equations are most commonly used in population modeling and are capable of generating fluctuations, including those that are complex in shape [44] . These are the Hutchinson, Nicholson, and Gopalsamy models, which have several modifications [45] . The most common Andronov-Hopf bifurcation reflects a mild loss of a stable regimen and describes the emergence of the stable cycle well in the vicinity of a stationary point that loses its stability [46] . Modifications of the Hutchinson model are used to simulate the harmonic fluctuations established in an isolated population without effects of external factors or overshoot, which is a situation where a rapidly growing invader population quickly exceeds the niche capacity K: t m < t c , N(0) < K : N(t m ) > K, N(t c ) < K. Wright [47] proposed Eq. (4), which is based on Hutchinson's ideas and is an extension of the Verhulst model with the lag τ in the regulation of reproduction: The decaying fluctuation mode = K it can be obtained with Eq. (4) when 1 < rτ < π/2, N(0) < K. However, Eq. (4) cannot describe a scenario where fluctuations stop with a decrease in population size to the minimal stable level = L, even by using the independent decrement. At rτ = π/2 + ε, the focus point loses its stability and a bifurcation occurs to give birth to a new cycle [48] . Harmonic fluctuations with a regular damping are rarely observed in ecological reality. Cycles known for small rodents of northern regions are a sequence of increasing annual generation sizes with an abrupt maximum at the end of a period. The population returns to the minimum level L after a maximum [49] . Sawtooth-like outbreaks of insect pests are individual series of sporadic Λ-shaped peaks that tend to decrease in amplitude, with the second peak being maximal in many cases [50] . A scenario known as sawtooth-like oscillations in insect pest dynamics includes population size outbreaks with different-sized decaying peaks; it differs from other variants of outbreak development [51] . An external effect is necessary in the scenario in order to distort equilibrium and to induce the next transition mode. Figure 8 shows a series of decaying outbreaks observed for Lymantria dispar in Bulgarian forests. The regular series of 2020 ended abruptly without repetition over a decade because of the activity of the entomopathogenic fungus Entomophaga maimaiga, which was a new larval parasite [52] . The capability of sensing and avoiding dangerous fungi is acquired with time by some insects, such as adapted termites, during their coevolution. A scenario with short outbreaks, which is characteristic of psyllids in southern Australia, is a biocybernetic phenomenon with a different regulatory mechanism, a threshold-dependent development of phases in the process [53] . A threshold population size sufficient to start an outbreak can be achieved from the stochastic oscillatory mode in Psyllidae. The existence of a threshold is associated with the activity of superparasitic wasps. A psyllid outbreak ends with a transition to a prolonged aperiodic mode. Threshold effects are important as a cause of resource collapses in the case of harvested species [54] . A problem known for the practical applications of model (4) is that minimums in the vicinity of zero arise with the increasing amplitude in a nonharmonically shaped cycle that forms at higher rτ values: = 0 + ε. Prolonged deep nearzero minimums disagree with the population dynamics observed in species that show eruptive phases in the rapid development of population size outbreaks. The area of affected forest regions reflects the intensity of a pest outbreak. Lack of newly dead forests in reports of the forest service means that the local pest population occurs in a state of a stably low density, but does not suggest its extinction. Thus, the Hutchinson model loses the possibility of population interpretation of fluctuations with the increasing rτ > 3π because various stochastic disturbances come to play a substantial role in a prolonged state with a low population size, even when the population does not experience the Allee effect. To eliminate the drawbacks of the cyclic solution of Eq. (4), Gopalsamy [55] proposed a model where the actual capacity of the ecological niche is no longer static: (5) Model (5) was modified to irregularly increase the nonlinearities in the numerator and denominator to provide that a new peak starts from the position min higher than the stable population size threshold L and Ψ(N k (t -h) normalizes the oscillation decay rate: The niche capacity for environment-damaging pests is not static in this model. Equations without explicit fixation of levels, such as the achievable optimal population size C or the limit capacity K of the ecological niche, provide alternative modifications of Eqs. (4) and (5) . The blowfly equation was proposed to describe the data (larva and adult counts) from experiments with a laboratory fly population at a high fluctuation amplitude: where δ is the independent mortality index. Model (7) allows an aperiodic fluctuation mode at higher τ ( Fig. 9 ). Competition between larvae and adults was implied in the relevant experiments [56] . The blowfly equation provides only a phenomenological description of relaxation fluctuations of total population size, taking advantage of the properties of the Ricker function. Models (4) and (7) were generalized to obtain model (8) for the age structure, but the resulting model did not find broad application: (8) In place of using the Ricker exponential regulation function f(x) = axe -bx to modify model (7), the nonunimodal function f(x) can be used in order to formalize intraspecific competition without using fixed population sizes as competition or niche saturation thresholds. Unimodal variants of f(x) were considered in [57] . Separate consideration can be given to the modification that is based on Bazykin's idea of an effect of the critically low population size L in the form of (N -L) and introduces a lag, as in the Hutchinson model: (9) Cyclic changes observed in reality for species abundance can be associated with the presence of reproductively isolated populations. Groups of this kind are impossible to distinguish in harvesting statistics, as was noted for the Volga sturgeon [58] . Overall statistics are misleading when obtained for a structurally complex population. To model the crisis effects that occur during rapid invasion, methods used to describe the delayed dynamics should be extended to include not only the intraspecific regulation, but also counteraction, which may also change in a threshold-dependent manner, for example, at a different initial infectious dose of influenza virus [15] . In the simulation methodology under study, the lag parameter in the dynamics of populations reproducing in a competitive environment should be understood as an aggregate characteristic that arises from a superposition of several biosystemic factors and particular concomitant processes. The temporal aftereffect as a component of regulatory functions of the models developed does not characterize a particular biological species, but reflects the systemic interactions that involve its surrounding and are often indirect rather than direct. Replication-and adaptation-related time indices were isolated separately for the lags associated with activation of different immunity components in model (3) of virus infection at various intensities of the immune response. The ontogenetic lag in the development of an adult is not the only temporal factor and is not equally significant in all situations. Lags that regulate reproductive activity are of greater importance. The parameter τ in model (2) of oscillation disruption is a result of the interaction that forms between the species and its particular biotic environment in a particular region. Contributions to the actual lag used in a population model can be made by ontogenetic properties, time periods of developmental stages, the restoration rate of vital resources, migration duration, etc. The time it takes for a parasite or a pathogen to develop a response is important in the case of invasion. The lag is an individual-dependent characteristic of all activation reactions in the case of immunity-virus interactions in the body, and its mean value may vary among local subpopulations. As was shown in [59] , the reproductive component rN(t) of a population model is more adequate without a lag, while the τ value is significant in the spatial component. Delayed equations of the N(t -τ) form are broadly used not only to describe the fluctuation dynamics of a population [60] or to simulate the immune response (modern developments of studies by Marchuk et al. [61] ), but also to study the intracellular processes, such as protein synthesis, nucleic acid replication, etc. When alternative RNA splicing was simulated [62] , the lag was similarly understood as an aggregate parameter of two main components, the total time of pre-mRNA splicing to produce mRNA and the time of mRNA transport from the cell nucleus into the cytoplasm. When a strictly controlled process is to be simulated, it is important that trajectory behavior remains predictable in computational experiments because the solution may become sensitive to minor disturbances in the initial conditions (N(0) ± ε) and the choice of the background function φ 0 for the time range φ 0 (ζ), -τ ≤ ζ < 0. Based on the above principles of the role played by the regulation with an aftereffect, modifications of the proposed equations will include the internal regulation of reproductive activity with τ (possibly, using threshold states of species abundance) and, separately, the functional of delayed external counteraction by the environment. Let the logarithmic regulatory function ln(K/N) be used in the model as a basic function for further special modifications at a stationary niche Decaying and cyclic fluctuations can be obtained in model (10), but another behavior of the equation solution is interesting from the ecological viewpoint. This scenario, where balanced equilibrium is achieved via oscillations in the solution of Eq. (10), differs from the scenarios describes by known models and is possible when N(0) is low and τ is sufficiently high (Fig. 10) . The parameters used in a comparative computational experiment were selected so that coinciding maximums were obtained for the solutions of Hutchinson Eq. (4) and basic model (10) . In the case of rapid invasion, an initially small group of the invader species can quickly, within a limited short period of time after invasion, increase to a huge population size (N(t K ) @ K) that is one order of magnitude higher than the final value achieved in equilibrium with the environment. The rapid growth is followed by a similarly rapid decrease in population size and oscillations around the equilibrium state, and the resulting population size, which slightly oscillates within a certain range, does not exert pressure on the environment. The concentration dynamics observed for the bacteriophage in experiments after its introduction matches this scenario, while a qualitatively different scenario is demonstrated by antagonistic bacteria. An outbreak that follows scenario (10) is possible for an initial group whose N(0) is substantially lower than K. There are other experimental examples where an outbreak ended. When the delicatessen scallop Chlamys farreri was grown artificially in marine cages, its production increased rapidly (Fig. 11) . Virus infection caused a six-fold decrease in production in 1997 and no measure could restore the production to its level of 1996 [63] . Parasitic wasps have been repeatedly introduced as a means to control pests, but all attempts at establishing their stable populations failed. A decrease in abundance was recently observed even in invasive mollusks of the family Dreissenidae, which occupied North American water bodies and became an important economic problem in the United States. Competition arose between the invasive species Dreissena polymorpha and D. rostriformis bugensis and was so intense that their total biomass density is lower than the density in a niche occupied by one of the species [64] . A single outbreak in the abundance of an invasive species is described by the model. Repeated abundance peaks [65] that arise after the first outbreak are characteristic of insect forest pests. A transition to sawtooth-like fluctuations is observed in pests, as mentioned above. The scenario with short Λ-shaped outbreaks in aphids and psyllids requires other attractor modifications [66] . The dynamics of model (10) do not describe the case where acute virus infection develops into chronic disease because a high virus load is maintained for a long period of time in the vicinity of unstable equilibrium (a developmental threshold) in this scenario. Pressure exerted by the environment may linearly depend on the previous state of the invader. Let us allow for a delayed response of the biotic environment in the modification of model (10) with a counteraction factor, but sacrifice the convenient property , N(0) > 0, N(t) ≥ 0: (11) where δ is the factor of mortality due to parasites or virus infections. The behavior of model (11) in a computational experiment shows extinction of the population only after a repeated outbreak. The second maximum exceeds the first peak in the solution (Fig. 12) , but the second minimum is outside of the allowable range of values, and when model computations are performed at N(t > τ) < -ε the computational experiment is completed by the instrumental algorithm with a fixation of extinction of the population. Model (11) phenomenologically reproduces the dynamics observed in Gause's experiment with ciliates. However, this is an artificial laboratory scenario observed in vitro, where the death of the invading predator is determined by rapid exhaustion of food resources. The time point t 0 of extinction N(t 0 ) = 0 depends on N(0) in model (11) , and only a single max- Argentine ant Linepithema humile is the best known undesirable invasive species whose populations were described many times to become extinct after an outbreak. The aggressive invader formed large colonies in New Zealand, but then the colonies collapsed suddenly, and autochtonous ant species restored their ranges after being displaced by the invader [67] . The scenario where invader populations intensely exhaust their necessary resources, encounter a contraction specified linearly with a lag, and eventually become extinct is important in the context of artificial control of undesirable invaders. Certain species (e.g., parasitic hymenopterans) are now often used to control other invaders biologically [68] , but complete control is rarely achieved. It is interesting to consider the South Korean experience of long-term fight against the American mulberry moth Hyphantria cunea, which is a quarantine pest [69] . In this scenario, the counteraction exerted by the environment (or an introduced species) may be determined by the state of the invader population at its peak size and at the time point when the population size has already passed its maximum as a result of independent factors. External pressure may be more complex in its organization than the proportional elimination -δN(t -τ). The response of introduced parasitic wasps is determined by the swarming of their prey and pathogenic viruses reproduce faster in dense agglomerations. Epizootics often affect intensely reproducing invaders as * N their population size approaches a certain threshold J < K, but long before the dangerous invader becomes capable of exerting its harmful effect on the environment. The third modification of the model addresses the most interesting situation, where a delayed counteraction to a rapidly growing population is nonlinear and depends on a threshold: (12) where J is the threshold at which the aggression of the biotic environment increases and m is the strength of the environmental response to the invader density. A regulatory mechanism that is complicated in terms of biological cybernetics unexpectedly yields a more flexible evolutionary scenario. The population size increases from the initial small value N(0) < J according to the common logarithmic law at the first stage of the process. However, in place of N(t) → K, the population size reaches a transition oscillatory mode far sooner than in model (4) or (7) . The growth shows neither gradual deceleration after a point of inflection nor a stop at a stage with the overshoot N(t K ) > K in new model (12) ; however, a crisis starts abruptly at t > 2τ (Fig. 13) . In a generalized sense, a model with the constant delay factors τ and h is a balance of the functions Φ, which describes the regulation of reproduction efficiency, and Ψ, which describes the counteraction of the biotic environment, with an important role being played by the ratio of m and d: When the threshold is approached (N(t) → J), the mortality rate greatly exceeds the birth rate and the population passes a bottleneck. As a result of passing the crisis, the solution trajectory tends to equilibrium at a level lower than J; the parameter m determines the crisis severity. Oscillatory dynamics are preserved at q = 0. The model is applicable to the case of small invader groups; the crisis becomes critical V-shaped when the threshold J is approached; the solution of Eq. (12) shows a form of logarithmic growth at N(0) > J. A substantial τ value does not lead to a stable cycle in the model. The threshold J in the model is a mirror reflection of L included in model (9) as the lower limit of population size from the model that Bazykin designed to describe the Allee effect. The fac- tors responsible for the origin of threshold population states and barriers in ecodynamics were considered in [70] . An example where a threshold arises from above and is due to activity of a superparasite was provided in [71] , which focused on the mechanisms that regulate reproduction in eucalyptus psyllids. A threshold that the population passes in a discrete model with a bifurcation following the outbreak does not arise in the pest, but arose in its natural enemies due to attacks by superparasites [72] . The specific parasitic efficiency of introduced wasps of the family Trichogrammatidae was found to increase with the increasing density of moth eggs [73] . As a result, a dangerous outbreak can occur in an invasive population with a high reproductive potential only when the condition N(t) > J is met, but intense counteraction does not lead to a critical minimum in population size in this case. The population successfully passes the bottleneck because the pressure that the environment exerts on the population and that depends on its previous state weakens when the current population size becomes substantially lower than the threshold. The model with complicated regulation of environmental counteraction provides an example of adaptation dynamics associated with passing through a deep crisis, as in the situation where a bacteriophage is introduced in a colony of bacteria that possess the CRISPR antivirus mechanism. The scenario obtained in the new model is applicable not only to parasitehost coevolution. Insect pests similarly adapt to resistant cultivars of genetically modified plants [74] . Bacteria and myxomycetes fight a chemical war, which led to the discovery of penicillin, but early antibiotics lost their efficacy with time. Invasive species often distort the balance in a trophic chain [75] to trigger a domino effect and to cause dramatic fluctuations in other species. The above equations can be used as components of systems designed to model direct trophic interactions. The explicit effect P of an antagonistic species can be introduced in Eq. (12) in place of q: (13) To describe the refugium effect with the possibility for a small relic population to get extinct, a power of one-third can be used for the effect of the lower threshold L: (14) The dynamics P(t) of a competitive interacting species is interesting to describe using fluctuations with another regulatory mechanism, as will be discussed in a forthcoming article. Models (12) and (14) [76] . Cases where excess nonlinearity is better to avoid in biosystem models employed in computational analyses were described in [77] because this leads to chaotization following the Feigelbaum scenario with an infinite cascade of bifurcations with the period p = 2 i , i → ∞ of the cycle ϕ(x*) n = ϕ(x*) n + p in the iterations x n + 1 = ϕ(x n ). Iteration models may include several excess effects, such as the disruption of unusual attractors, the formation of fractal boundaries in attraction domains, and the appearance of periodicity windows with chaos-cycle transitions and cycles following the order of Sharkovskii's theorem [78] , rather than ecological principles. An essential biological interpretation is necessary for each mathematical effect that arises in the proposed model. Hybrid models with event-based time representation are expedient for situations of a controlled anthropogenic effect or a crisis due to bioresource exploitation [79] because the control logic must be formalized. CONCLUSIONS Models (10), (11) , (12) , and (14) were proposed and used to consider the abrupt changes in population size in special cases of population processes. Such changes were sometimes classed with eruptive ecological processes or described as boom and bust dynamics in modern terms [80] . Example ecological situations that best match the solution behavior were described for each of the models. Gradient and eruptive processes are recognized in biological communities according to a published classification [81] . Gradient changes, such as succession, are easy to explain. A case where the situation changed in favor of one of the species in a manner independent of that species can be reflected by an increase in K with a subsequent monotonic, though rapid, increase in abundance. Eruptive and nonmonotonic dynamics are actually diverse, rather than being limited to a cyclic alternation of growth and decline phases of competing species. Such cases present a separate problem in methodology of biosystem modeling. The study [82] considered 16 situations where a rapid spreading of an invasive species was followed by a similarly rapid collapse. A collapse was understood as a situation where the population size decreases by 90% of its optimum and remains low for at least 3 years. Models (10) and (12) were used to consider the variants where rapid growth phase alternate with deep crises, up to risk of sudden extinction in model (14) . Using data from Gause's experiments and experiments with bacteriophages, delayed equations were shown to be efficient to use not only for a classical problem of modeling long-term stable endogenous fluctuations in population size [83] or decaying sawtooth-like outbreaks in moths. With the threshold regulatory functional Ψ[N(t -h)] included in the models with N(t -τ), the properties of extreme forms of invasion dynamics and scenarios with rapid adaptation can be described. Delays in regulation are important to analyze in various processes [84] that have several distinct developmental stages and conditions for a transition between stages. The previous state-dependent linear external effect Ψ = -δN(t -τ) can lead to extinction of the population in model (11) . In modification (12), a rapidly growing population enters a deep crisis before achieving the size that is theoretically possible in its ecological niche. However, the regulation is more complex in this scenario and this addition makes it possible to successfully pass a bottleneck stage, as E. coli colonies demonstrated in experiments. The models developed in this work describe the scenarios where dynamic modes change without bifurcations; such transitions in scenarios do not require the parameter values to be adjusted in the course of a model experiment. To simulate a process, it is important to understand which ones of the observed changes can be modeled using bifurcations of attractors, that is, the formation of topologically inequivalent phase portraits of the trajectory in response to disturbances of an internal parameter. Parameters of nonlinear models differ in effect on the phase portrait topology. Shifting the vector of equation parameters is the simplest method, but is not always justified biologically. Internal characteristics of biosystems evolve in a concerted manner. A stepwise change in characteristics seems unrealistic in the majority of cases (with the exception of rapidly arising mutations in viruses and prokaryotes), and different mathematical tools are consequently developed in my series of works. The complex threshold-dependent regulation that was considered using model (12) and is due to the fact that the regulator species has a limiting factor is not a unique phenomenon. It is quite possible to encounter a real situation where one species disrupts the threshold (or creates a new one) without being involved in direct interactions. Bacterial pneumonia may develop because of virus infection when the population of alveolar macrophages is exhausted [85] . Multilevel regulatory schemes that act to control reproduction in species with a high reproductive potential are interesting for biological cybernetics, but it is difficult to provide direct mathematical descriptions for their components. The reproduction rate and the life cycle duration of a species may be found to belong to different or even incomparable time scales. A simplification to a closed predator-prey system leads to N = 0 based on Gause's experiments. Evolutionary adaptation may create unexpected barriers that prevent the total dominance of an efficient predator or parasite [86] . As an example, the subvirus parasites virophages utilize large amoeba-infecting viruses for their reproduction and thus increase the viability of amoeba [87] . Bacteria of the genus Bdellovibrio are parasites of other bacteria and may find application in medicine in the future [88] . One of the most complex systems is known for tropic molds [89] and includes tens of species, which compete and chemically suppress each other. Sharp competition for limited resources acts as a factor of intense speciation to generate many endemic species of one family not only in microorganisms, but also in other organisms, as is the case with 35 fish species of the family Gobiidae in the Caspian Sea [90] . Many discoveries of applied significance can be expected from studies of the evolutionary struggles in the microbial world. As an example, archaea seem to acquire a separate antiviral mechanism independent of the bacteria CRISPR/Cas system [91] . Types of lags should be isolated when developing systems of equations according to the above principles. Reproductive lags, regulatory lags as a factor of a delayed response of the environment, and adaptationrelated lags can be recognized. It takes time for a population (or an organism) to develop adaptive responses, and this lag may change during an invasion. The effect exerted by the biotic environment is dynamically adjusted in the most interesting invasion cases. The population size is not the only parameter that is affected, as was observed when North Atlantic cod stocks collapsed because of a systemic, though minor, error in estimating the total allowable catch of adult cod [92] . The function regulating reproduction efficiency undergoes transformation and the parameter m in model (12) may depend on t. There are processes where the three lag factors compete, e.g., different mechanisms of the immune response differ in activation and inhibition times, as is also observed for virus replication cycles. Dysregulation and especially hyperactivation of the immune response increase the pathogenesis of COVID and HIV infection [93] . Cyclic solutions are not important in modeling the fight between a virus and immunity. We note that logic-algorithmic methods make it possible to include even a greater number of eventdriven change factors in models of confrontation between spatially nonuniform populations [94] . Higher-density fronts often form when an invasive species spreads through a new area, as is also observed in common microbial communities [95] . Lag factors can be included in cellular automaton algorithms [96] , but with special rules [97] . Software implementations of intellectual cellular automata with a variable cell density in a grid and a stochastic component in the algorithm are a promising field of current research and can be used to predict the spatially nonuniform processes that take place in biosystems. The Struggle for Existence Stochasticity of Dynamic Systems Il'ichev Delay Differential Equations and Applications Stability and Complexity of Ecosystems The Hopf Bifurcation Stability and Oscillations in Delay Differential Equations of Population Dynamics Mathematial Problems in Cybernetics Encyclopedia of Entomology The author declares that he has no conflict of interest. This work does not contain any studies involving animals or human subjects performed by any of the authors.