3.0cm3.0cm2.5cm1.5cm

# Modelling Ethnogenesis

Alexey Piunovskiy
Bakhti Vasiev

Department of Mathematical Sciences, University of Liverpool, L69 7ZL, UK.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Following the ideas of L.N.Gumilev, we introduce the mathematical model of
ethnogenesis which describes the dynamics of subgroups in the developing
polity in terms of ordinary differential equations. The bust dynamics
associated with the rise and fall of civilisations is modelled as an
excitation process, which is the non-linear phenomenon, well known in
mathematical biology. We consider deterministic as well as the stochastic
version of the model. We also expand the model to study the interaction
between two polities undergoing ethnogenesis. Investigation is performed using
analytical methods as well as numerical integration (i.e. MATLAB simulation).

|  Keywords: Ethnogenesis, Population Dynamics, Dynamical System, Excitation
---|---
|  AMS 2020 subject classification: | 37N25, 37M05, 92B05, 92D15, 92D25

##  1 Introduction

History of mankind in the last 5 thousand years can be viewed in terms of
rises and falls of civilizations. L.N.Gumilev considered each civilization as
a manifestation of particular ethnos, which, under certain conditions, appears
and builds civilization, but in course of time gets old and dies, causing for
the associated civilization to disappear. According to Gumilev’s theory [3, 4]
the driving force for formation of a new ethnos comes from a certain group of
people, whom he calls ”passionaries” (the term which is also used in [10]) or
people with drive. This group is considered as a fraction of population
habitating certain geographic territory, who express high level of passion and
lead their fellows (i.e. tribesmen) forming the rest of the population to
expand and to build new society (civilization).

According to Gumilev [3, 4] the formation and death of civilization can be
described by the dynamics of civilisation’s ”passionary tension” or drive,
which he illustrated by a ”bust” curve shown in Figure 1. It starts with a
growing phase (rise of civilization) followed by plateau (ackmatic phase),
fast decline (breaking phase), slow decline (inertial obscuration) and low
level tail (obscuration or regeneration-relict). The entire process, which
according to Gumilev takes about 15 centuries, can be considered as a response
to a disturbance, caused by initiation of a small fraction of passionaries.
Such response is known in physiology as excitation [7, pp.239-242], that is
when a small perturbation to the system results in a full-sized response.

Figure 1: Graphical illustration of the evolution of an ethnos (civilization)
[3, p.240].

The aim of this study is to develop a mathematical model which would explain
the bust dynamics (which is evident from Figure 1) exhibited in course of
ethnogenesis. We note that mathematical study of ethnogenesis was also aimed
by other researchers. However, the mathematical models developed in their
works (see, for example, [5, 10], do not appropriately reproduce the dynamics
of ethnogenesis, as they didn’t consider the association of the bust dynamics
with the type of system’s nonlinearity known as excitability. The excitable
dynamics in the model describing the ethnogenesis is the main theme of the
current work. We will build a mathematical model of ethnogenesis, which is
based on the statements underlying Gumilev’s theory.

  * •

According to Gumilev there are three main subgroups in the population having
different behavioural patterns and affecting the evolution of the ethnos. The
driving force for the growth of populations and further improvements resulting
in the rise of a civilization comes from passionaries, whose idealistic
motivations are grounded on altruism. The bulk of population is represented by
harmonious individuals who work on preservation of the current state of the
ethnos. There is also a destructive group of individuals, called
subpassionaries (or people with negative drive – vargants, soldier tramps,
degenerates), who are as active as passionaries but whose actions are based on
egoism rather than altruism.

  * •

The formation of a new ethnos (or civilisation) is associated with initiation
of passionaries. Gumilev describes conditions under which passionaries appear
and take over the population, but we will not go into details of these
conditions and will postulate that at a certain time, a small fraction of
population is already represented by passionaries.

We will design a few models describing the ethnogenesis. In the first model we
will consider the population as consisting of two subgroups, namely, the
passionaries and the remainder of the population. This two-variable model will
let us to identify possible interactions between these subgroups which allow
the bust dynamics in course of ethnogenesis. In the second version of the
model we will consider all three subgroups and analyse the interactions
between them which is consistent with the observed dynamics during
ethnogenesis. In the follow up steps of our research we will use the three-
variable model for the study of the impact of noisy environment to the
ethnogenesis. Particularly we show that the noise is amplified by
nonlinearities in the model. Finally, we extend the three variable model (with
noise) to consider interaction between two ethnogenetic processes taking place
simultaneously, with a certain time lag. This allows us to model conflicting
civilisations and to identify conditions when one of them takes over the
other.

##  2 Two-Variable Model

As a starting point we will consider the two-variable model:

$\left\\{\begin{array}[]{l}\dot{x}=xf(x,y)\\\
\dot{y}=yg(x,y)\end{array}\right.$ (1)

where variable $x$ represents the size of the subpopulation formed by
passionaries while $y$ is the size of the remaining population (which includes
both harmonious people and subpassionaries). Equation (1) is commonly used for
modelling population dynamics in biology, where $x$ and $y$ are considered as
the sizes of two biological species. Using linear approximation of functions
$f(x,y)$ and $g(x,y)$ we get:

$\left\\{\begin{array}[]{l}\dot{x}=x(a_{0}+a_{1}x+c_{1}y);\\\
\dot{y}=y(b_{0}+b_{1}y+d_{1}x),\end{array}\right.$ (2)

which is a generalised representation of the Lotka-Volterra model [11].
Commonly this model is considered under the following condition for model
parameters: $a_{0}>0$, $b_{0}>0$ while $a_{1}<0$ and $b_{1}<0$ preventing
unlimited growth of populations. After nondimensionalisation this equation is
commonly transferred into

$\left\\{\begin{array}[]{l}\dot{x}=x(1-x+\beta_{1}y);\\\ \dot{y}=\gamma
y(1-y+\beta_{2}x),\end{array}\right.$ (3)

where $\gamma$ defines the relative rate of change of $y$ with respect to $x$
[7, p.119]. If $\beta_{1,2}=0$ the two equations are detached and the both
species exhibit the logistic growth. Furthermore, depending on the signs of
these two parameters the model reproduces three types of interactions between
populations $x$ and $y$, namely, predator-pray ($\beta_{1}>0$, $\beta_{2}<0$),
symbioses ($\beta_{1}>0$, $\beta_{2}>0$) and competition ($\beta_{1}<0$,
$\beta_{2}<0$) [11]. The model allows four equilibria and in case of
competitive Lotka-Volterra model ($\beta_{1}<0$, $\beta_{2}<0$) all four
equilibria are meaningful and correspond to non-negative sizes of populations.
These are the trivial equilibrium $(x=0,~{}y=0)$, the extinction of
$y$-population ($x\neq 0$, $y=0$), the extinction of $x$-population ($x=0$,
$y\neq 0$) and the co-existence ($x\neq 0$, $y\neq 0$). Depending on model
parameters, the solution of the system converges either to the co-existence of
species or to the case when one of them becomes extinct [7, pp.104-126].

In our case we would like for one of the variables, say $x$ to represent the
size of the subpopulation formed by passionaries, while $y$ is the size of the
remaining part of population (which would include harmonious people and
subpassionaries). As we expect to observe the excitation dynamics, and the
excitation in biology is known to be a non-linear process, we will need more
than just a linear expansion of functions $f(x,y)$ and $g(x,y)$ in (1). The
simplest way is to add one quadratic term into the first equation which
transforms our system to

$\left\\{\begin{array}[]{l}\dot{x}=x(a_{0}+a_{1}x+a_{2}x^{2}+c_{1}y);\\\
\dot{y}=y(b_{0}+b_{1}y+d_{1}x).\end{array}\right.$ (4)

There are up to six equilibria in this system of which two are always real:

$(x_{1},y_{1})=(0,0),(x_{2},y_{2})=(0,-b_{0}/b_{1}).$

These two equilibria are located on the vertical axis. Two more equilibria (if
real) are located on the horizontal axis:

$(x_{3,4},y_{3,4})=\left(\frac{-a_{1}\pm\sqrt{a_{1}^{2}-4a_{0}a_{2}}}{2a_{2}},0\right).$
(5)

And, finally, two remaining equilibria (if real) are represented by the points
of intersection of the line $b_{0}+b_{1}y+d_{1}x=0$ with parabola
$a_{0}+a_{1}x+a_{2}x^{2}+c_{1}y=0$. To set an excitable kinetics in the system
(4) we make sure that the equilibrium $(x_{2},y_{2})$ is stable and located in
the vicinity of the parabola $a_{0}+a_{1}x+a_{2}x^{2}+c_{1}y=0$. The excitable
dynamics becomes more evident after nondimensionalisation of the system (4)
(so that $a_{2}=-1$) and transferring it into the following form:

$\left\\{\begin{array}[]{l}\dot{x}=x((1-x)(x-\alpha)+\beta_{1}(y-y_{0}));\\\
\dot{y}=\gamma y(y_{0}-y+\beta_{2}x),\end{array}\right.$ (6)

where parameters $\gamma$, $\beta_{1}$ and $\beta_{2}$ have the same meaning
as those in the equation (3), new parameter $\alpha$ defines the excitation
threshold of the system and $y_{0}$ defines the location of the equilibrium
$(x_{2},y_{2})=(0,y_{0})$ which is stable if $\beta_{1}$ and $\beta_{2}$ have
opposite signs. Note, that the parameters $\alpha$, $y_{0}$ and $\gamma$
should all be positive. It looks that, in any case, the system (6) with
positive initial conditions admits a unique bounded solution on the infinite
horizon $[0,\infty)$. For $\beta_{1}\leq 0$, $\beta_{2}\geq 0$, this will
follow from the investigation of the three-variable model (7). Concerning the
equilibria in the system (6), we note that the trivial steady state
$(x_{1},y_{1})=(0,0)$ is unstable (saddle). Furthermore, if the parabola
$(1-x)(x-\alpha)+\beta_{1}(y-y_{0})$ has real and non-negative roots, then the
system (6) has two meaningful equilibria located on the horizontal axis (same
as given by (5)) and for the concave up parabola the equilibrium which is
closer to the origin (smaller $x$-coordinate) is the unstable node, while the
other one is a saddle. For simplicity, we consider the cases when the
nullcline represented by the parabola $(1-x)(x-\alpha)+\beta_{1}(y-y_{0})=0$
doesn’t intersect the one given by the line $y_{0}-y+\beta_{2}x=0$ and
therefore we don’t have any extra equilibria.

Figure 2: Illustration of the excitable dynamics observed in the system (6).
Panel A: nullclines $\dot{x}=0$ and $\dot{y}=0$ are shown in blue and a set of
phase trajectories with starting points $x=0.05,0.1,...,0.5$ and $y=y_{0}$ are
shown in red. Panel B: time evolution of both variables, $x$ (blue line) and
$y$ (red line) from the initial condition $x=0.1$, $y=y_{0}$. Parameters
values: $\alpha=0.02$, $y_{0}=0.05$, $\beta_{1}=-1/3$, $\beta_{2}=2.5$,
$\gamma=0.1$.

The dynamics in the system (6) is illustrated by Figure 2. Null-clines of the
system are shown in blue on panel A. These null-clines indicate the excitable
nature of the system (6) and this is illustrated by a set of phase
trajectories (shown in red) coursed by the perturbation of the system from its
stable equilibrium $(x_{2},y_{2})=(0,y_{0})$. Any perturbation from this state
results in the relaxation of the system back to this equilibrium. However if
the perturbation is above certain threshold, for example, if the initial value
of $y$ is equal to $y_{0}$ and the initial value of $x$ is above $\alpha$,
then the perturbation increases further before the system relaxes back to the
stable equilibrium $(x_{2},y_{2})$.

The dynamics of variables $x$ and $y$ over time for one of the phase
trajectories (starting from the point $x_{0}=0.1$ and $y_{0}$) is shown in
panel B. Here we see that the both variables increase and then decrease over
time. We note a relatively fast dynamics of variable $x$ (passionaries) with
the duration of the spike being about 60 time unites. We also note a slow
relaxation of variable $y$ (the rest of population) which gets back to its
equilibrium value with relaxation time of about 100 time units.

In order to scale the model time units to real time we take into account that
according to Gumilev [3] the duration of passionary spike is about 900 years
which should correspond to 60 time units in the model. Thus one model time
unit corresponds to 15 years. Furthermore, according to Gumilev [3] the
fraction of passionaries can only be up to 5-7% of the entire population. For
the set of parameters values used to produce Fig. 2 both variables have values
roughly in the range (0,1). If we consider variable $x$ as representing 1% of
$y$, the total number of nonpassionaries, then at the top of the spike (which
takes place at $t\approx 15$ and where $x\approx 1$ and $y\approx 0.2$) the
passionaries constitute about 5% of the entire population. In other words,
$x=1$ corresponds, e.g., to the absolute value of $K=10,000$ passionaries,
while $y=1$ corresponds to $100K=1,000,000$ nonpassionaries. One can certainly
take other values of $K$. Note also that on the time interval $(20,30)$ of the
most rapid growth of nonpassionaries (constituting the main part of the
population) their amount doubles. That corresponds to doubling time
$\frac{10\cdot 900}{60}=150$ years which is roughly in line with observations
on the maximal growth rate of human populations.

Figure 3: Illustration of the excitable dynamics in the system (6). Panel A:
nullclines $\dot{x}=0$ and $\dot{y}=0$ are shown in blue and a set of phase
trajectories with starting points $x_{0}=0.05,0.1,...,0.5$ and $y_{0}$ are
shown in red. Panel B: time evolution of both variables, $x$ (blue line) and
$y$ (red line) from the initial condition $x=0.1$, $y=y_{0}$. Parameters
values: $\alpha=0.02$, $y_{0}=1$, $\beta_{1}=1/3$, $\beta_{2}=-2.5$,
$\gamma=0.1$.

For the dynamics illustrated by Figure 2 it is essential that $\beta_{1}<0$
and $\beta_{2}>0$, that is, passionaries are suppressed by the rest of the
population while their own impact to non-passionaries is positive. It appears
that such relationships between passionaries and non-passionaries is not the
only one allowing excitable dynamics. An alternative case when $\beta_{1}>0$
and $\beta_{2}<0$, that is, passionaries are activated by non-passionaries
which are in turn suppressed by the passionaries can also result in the
exciatable dynamics. This scenario is illustrated in Figure 3 which similarly
to the Figure 2 shows the dynamics of the system (6) but with swapped signs of
the parameters $\beta_{1}$ and $\beta_{2}$ and with $y_{0}=1$ (rather than
$y_{0}=0.05$ in Figure 2). Panel A shows null-clines of the system and typical
phase trajectories obtained from the over-threshold perturbation of the steady
equilibrium $(x_{2},y_{2})=(0,y_{0})$. Perturbation is made by an increase of
the $x$-value over the threshold, $\alpha=0.02$. Time dependence of the
variables $x$ and $y$ for one of the phase trajectories is shown on panel B.
Now we see that $x$ and $y$ change in the opposite directions: initially $x$
is increasing (for $t<28$) and $y$ is decreasing (for $t<44$) and later the
both variables inverse their rate of change. Furthermore, here $y_{0}=1$ and
the entire dynamics can be seen as change in fraction of passionaries over
time for the population of a roughly constant size. Similarly to the case
shown in Fig. 2, model time unit corresponds to 15 years.

There is an important difference between the dynamics shown in Figures 2 and
3. In Figure 2 the spike in variable $x$ is followed by the spike in the
variable $y$, that is, an increase in number of passionaries is followed by
the increase of the size of remaining population. As for the dynamics shown in
Figure 3 we notice that the increase in number of passionaries is followed by
the decrease of the size of remaining population. We know that an increase in
number of pasionaries results in the expansion of the polity and
correspondingly to the growth of the population. Hence the dynamics shown in
Figure 2 looks natural if the variable $y$ represents harmonious people.
However if the variable $y$ is associated with subpassionaries then the
dynamics in Figure 3 is not impossible as subpassionaries may be suppressed by
the passionaries. Up to now the variable $y$ was considered as including both,
harmonious people and subpassionaries. In order to consider their dynamics
separately we will modify our model by allocating variables to each of these
two subpopulations.

##  3 Three-Variable Model

Gumilev in his theory of ethnogenesis considers three types of individuals who
constitute ethnos and whose behaviour has an impact to the ethnogenetic
process. To follow this concept we extend the two-variable model described by
the system (6) by including extra variable $z$, so that the variables, $x$,
$y$ and $z$, represent the sizes of subpopulations of passionaries ($x$),
harmonious people ($y$) and subpassionaries ($z$). Furthermore, we will
presume that the dynamics of subpassionaries is similar to that of
passionaries that is, their rate of change has quadratic dependence on their
own sizes. However passionaries and subpassionaries differ by their
relationships with harmonious people and each others. So, our three-variable
system can be represented as the following:

$\left\\{\begin{array}[]{l}\dot{x}=\gamma_{1}x[(1-x)(x-\alpha_{1})+\beta_{12}(y-y_{0})+\beta_{13}(z-z_{0})];\\\
\dot{y}=\gamma_{2}y(y_{0}-y+\beta_{21}x+\beta_{23}(z-z_{0}));\\\
\dot{z}=\gamma_{3}z[(z_{0}-z)(z-\alpha_{2})+\beta_{31}x+\beta_{32}(y-y_{0})]\end{array}\right.$
(7)

with the initial condition $x(0),y(0),z(0)>0$. Since its right-hand part is
locally Lipschitz, the system (7) admits the unique local solution [9, Theorem
2.2.]. It looks that in general it can be extended to the unique bounded
solution on the infinite horizon $[0,\infty)$.

Excitable dynamics can be observed in the system (7) under different kinds of
interactions between the variables. In general, the steady states of the
system (7) with $x=0$ are as follows:

  * •

$(x_{1},y_{1},z_{1})=(0,0,0)$;

  * •

$(x_{2},y_{2},z_{2})=(0,y_{0}-\beta_{23}z_{0},0)$;

  * •

$(x_{3,4},y_{3,4},z_{3,4})=(0,y_{3,4},z_{3,4})$, where $y_{3,4}$ and $z_{3,4}$
are the two solutions to equations

$\left\\{\begin{array}[]{l}y_{0}-y+\beta_{23}(z-z_{0})=0;\\\
(z_{0}-z)(z-\alpha_{2})+\beta_{32}(y-y_{0})=0.\end{array}\right.$

  * •

$(x_{5,6},y_{5,6},z_{5,6})=(0,0,z_{5,6})$, where $z_{5,6}$ are the two
solutions to equation

$(z_{0}-z)(z-\alpha_{2})-\beta_{32}y_{0}=0.$

Figure 4: Illustration of the excitable dynamics in the system (7). Panel A:
nullclines $\dot{x}=0$ and $\dot{y}=0$ at $z=z_{0}$ are shown in blue and a
set of phase trajectories with starting points $x(0)=0.04,0.07,0.1,0.4$ and
$y(0)=0.053$, $z(0)=0.05$ are shown in red. Panel B: time evolution of all
three variables, $x$ (blue line), $y$ (red line) and $z$ (black line) from the
initial condition $x(0)=0.07$, $y(0)=0.053$, $z(0)=0.05$. Parameters values:
$\alpha_{1}=0.03$, $\alpha_{2}=0.11$, $y_{0}=0.075$, $z_{0}=0.22$,
$\beta_{12}=-6$, $\beta_{13}=0.6$, $\beta_{21}=0.2$, $\beta_{23}=0.1$,
$\beta_{31}=0.5$, $\beta_{32}=0$, $\gamma_{1}=1$, $\gamma_{2}=0.7$,
$\gamma_{3}=0.2$.

In the case illustrated in Fig. 4, passionaries are suppressed by harmonious
people ($\beta_{12}<0$) and promoted by subpassionaries ($\beta_{13}>0$);
harmonious people are promoted by passionaries ($\beta_{21}>0$) as well as by
subpassionaries ($\beta_{23}>0$); subpassionaries are promoted by passionaries
($\beta_{31}>0$) and do not depend on harmonious people ($\beta_{32}=0$). For
the set of parameters values used in the simulation shown in Fig. 4, there are
eight real-valued steady states, with seven ones among them having non-
negative coordinates. All the steady states with non-zero component $x$ are
unstable. Among the six steady states enlisted above,
$(x_{1},y_{1},z_{1})=(0,0,0)$, $(x_{3},y_{3},z_{3})=(0,0.064,0.11)$,
$(x_{5},y_{5},z_{5})=(0,0,0.22)$ and $(x_{6},y_{6},z_{6})=(0,0,0.11)$ are
unstable, and the steady states $(x_{2},y_{2},z_{2})=(0,0.053,0)$ and
$(x_{4},y_{4},z_{4})=(0,0.075,0.22)$ are stable. In Fig. 4 the stable steady
states are shown with big blobs, and the unstable steady state
$(x_{3},y_{3},z_{3})=(0,0.064,0.11)$ between them is indicated as the short
line.

The excitation appears, starting in the neighbourhood of the stable point
$(x_{2},y_{2},z_{2})$$=(0,0.053,0)$ which is shown as the lower blob on the
vertical axis. We assigned the initial values $z(0)=0.05$ (to give a push from
the ‘cemetery’ $z=0$), $y(0)=0.053$; $x(0)$ varies from $0.04$ to $0.4$. If
the initial push $x(0)$ is below or slightly above the threshold $\alpha_{1}$,
the system quickly returns back to the state $(x_{2},y_{2},z_{2})$. But larger
(still small enough) initial perturbation results in the excitation leading to
the second stable point $(x_{4},y_{4},z_{4})=(0,0.075,0.22)$ shown as the
upper blob on the vertical axis. Depending on the value of $x(0)$, the
trajectory approaches the limit $(x_{4},y_{4},z_{4})$ either from below or
from above. Even a small over-threshold perturbation $x(0)=0.07>\alpha_{1}$
grows up to around $x=0.5$ before it relaxes back to
$\lim_{t\to\infty}x(t)=0$. Qualitatively, the picture is similar to that
presented in Fig. 2. Again, one model time unit corresponds to 15 years.

Figure 5: Illustration of the excitable dynamics in the system (7). Panel A:
nullclines $\dot{x}=0$ and $\dot{y}=0$ at $z=z_{0}$ are shown in blue and a
set of phase trajectories with starting points $x(0)=0.04,0.07,0.1,0.4$ and
$y(0)=0.075$, $z(0)=0.05$ are shown in red. Panel B: time evolution of all
three variables, $x$ (blue line), $y$ (red line) and $z$ (black line) from the
initial condition $x(0)=0.07$, $y(0)=0.075$, $z(0)=0.05$. Parameters values:
$\alpha_{1}=0.03$, $\alpha_{2}=0.11$, $y_{0}=0.075$, $z_{0}=0.$,
$\beta_{12}=-6$, $\beta_{13}=0.6$, $\beta_{21}=0.2$, $\beta_{23}=0.1$,
$\beta_{31}=0.5$, $\beta_{32}=0$, $\gamma_{1}=1$, $\gamma_{2}=0.7$,
$\gamma_{3}=0.2$.

As a special case, one can put $z_{0}=0$ leaving the other parameters the
same. For this set of parameters values, used in the simulation shown in Fig.
5, there is only one stable steady state $(x,y,z)=(0,0.075,0)$: all three
previous equilibria $(x_{2,3,4},y_{2,3,4},z_{2,3,4})$ now coincide. The system
again exhibits excitable kinetics: over-threshold perturbation of $x$
($x>\alpha_{1}$) grows up to around $x=0.2$ before it relaxes back to the
equilibrium.

An alternative dynamics for the system (7) is shown in Fig 6. Here the
interactions between the variables is slightly different from those used for
the dynamics illustrated in Fig 4. The difference is that harmonious people
are suppressed (rather than promoted) by subpassionaries ($\beta_{23}<0$) and
subpassionaries are promoted (rather than suppressed) by passionaries
($\beta_{31}<0$). For the set of parameters values used in the simulation
shown in Fig. 6 there are two stable steady states
$(x_{2},y_{2},z_{2})=(0,0.12,0)$ and $(x_{4},y_{4},z_{4})=(0,y_{0},z_{0})$:
the enumeration is in accordance with the expressions below equation (7). The
system at the steady state $(x_{4},y_{4},z_{4})$ exhibits excitable kinetics:
over-threshold perturbation of $x$ ($x>\alpha_{1}$) grows up to around $x=0.7$
before it relaxes back to the equilibrium $(x_{4},y_{4},z_{4})$. As we have
two stable states, the relaxation can bring the system to another steady state
$(x_{2},y_{2},z_{2})$, and this is observed in the system with slightly
different set of parameters values. Qualitatively, the dynamics illustrated in
Fig 6 is similar to that presented in Fig. 3.

Figure 6: Illustration of the excitable dynamics in the system (7). Panel A:
nullclines $\dot{x}=0$ and $\dot{z}=0$ at $y=y_{0}$ are shown in blue and a
set of phase trajectories with starting points $x(0)=0.05,0.1,0.2,0.5$ and
$y(0)=y_{0},z(0)=z_{0}$ are shown in red. Panel B: time evolution of all three
variables, $x$ (blue line), $y$ (red line) and $z$ (black line) from the
initial condition $x(0)=0.1$, $y(0)=y_{0}$, $z(0)=z_{0}$. Parameters values:
$\alpha_{1}=0.03$, $\alpha_{2}=0.1$, $y_{0}=0.075$, $z_{0}=0.6$,
$\beta_{12}=-0.06$, $\beta_{13}=0.6$, $\beta_{21}=1.25$, $\beta_{23}=-0.075$,
$\beta_{31}=-0.5$, $\beta_{32}=0$, $\gamma_{1}=2$, $\gamma_{2}=20$,
$\gamma_{3}=0.6$.

One can see that in Fig. 6 the number of subpassionaries dramatically
decreases before going back to the equilibrium. This perhaps does not often
occur in reality. The dynamic in Fig. 5 looks more reasonable. But the growth
of the passionaries subpopulation (from $0.07$ to $0.2$) is not as impressive
as in Fig. 4 (from $0.07$ to $0.49$). Therefore, in the further simulations,
we take the parameters values from the latter case (Fig. 4). If
$\beta_{12}\leq 0$, $\beta_{32}\leq 0$ the existence of the unique bounded
solution to the system (7) will follow from the investigation of the
stochastic version of the model: see Lemma 4.1 and its proof, especially,
Remark 7.2.

##  4 Ethnogenesis in Noisy Environment

In order to study the impact of noise to the ethnogenetic process we will
modify the three-variable model by adding extra (stochastic) terms to the
system (7). To justify the modification which we are about to impose, let us
consider the following change of variables:

$v_{1}=\ln x,~{}~{}~{}~{}~{}v_{2}=\ln y,~{}~{}~{}~{}~{}v_{3}=\ln z.$ (8)

In these variables the equations (7) transform into

$\left\\{\begin{array}[]{l}\dot{v_{1}}=\gamma_{1}[(1-e^{v_{1}})(e^{v_{1}}-\alpha_{1})+\beta_{12}(e^{v_{2}}-y_{0})+\beta_{13}(e^{v_{3}}-z_{0})];\\\
\dot{v_{2}}=\gamma_{2}[y_{0}-e^{v_{2}}+\beta_{21}e^{v_{1}}+\beta_{23}(e^{v_{3}}-z_{0})];\\\
\dot{v_{3}}=\gamma_{3}[(z_{0}-e^{v_{3}})(e^{v_{3}}-\alpha_{2})+\beta_{31}e^{v_{1}}+\beta_{32}(e^{v_{2}}-y_{0})],\end{array}\right.$
(9)

with the initial conditions $v_{1}(0)=\ln x(0)$, $v_{2}(0)=\ln y(0)$,
$v_{3}(0)=z(0)$. Note that $v_{1},v_{2}$ and $v_{3}$ may be negative and the
initial conditions $x(0),y(0),z(0)>0$ are assumed to be fixed.

The natural way to define the stochastic version is to introduce stochastic
differential equations

$\left\\{\begin{array}[]{l}d{V_{1}}=\gamma_{1}[(1-e^{V_{1}})(e^{V_{1}}-\alpha_{1})+\beta_{12}(e^{V_{2}}-y_{0})+\beta_{13}(e^{V_{3}}-z_{0})]dt+\sigma_{1}dW_{1};\\\
d{V_{2}}=\gamma_{2}[y_{0}-e^{V_{2}}+\beta_{21}e^{V_{1}}+\beta_{23}(e^{V_{3}}-z_{0})]dt+\sigma_{2}dW_{2};\\\
d{V_{3}}=\gamma_{3}[(z_{0}-e^{V_{3}})(e^{V_{3}}-\alpha_{2})+\beta_{31}e^{V_{1}}+\beta_{32}(e^{V_{2}}-y_{0})]dt+\sigma_{3}dW_{3};\\\
V_{1}(0)=v_{1}(0),~{}~{}~{}~{}~{}V_{2}(0)=v_{2}(0),~{}~{}~{}~{}~{}V_{3}(0)=v_{3}(0),\end{array}\right.$
(10)

which we understand in the sense of the Ito stochastic calculus [8]. Here
$W_{1}$, $W_{2}$ and $W_{3}$ are mutually independent standard Brownian
motions on the complete filtered probability space $(\Omega,{\cal F},({\cal
F}_{t})_{t\geq 0},P)$, and $\sigma_{1},\sigma_{2},\sigma_{3}>0$. After that,

$X=e^{V_{1}},~{}~{}~{}~{}~{}Y=e^{V_{2}},~{}~{}~{}~{}~{}Z=e^{V_{3}}$

will be the random processes representing the sizes of the subpopulations of
passionaries, harmonious people and subpassionaries respectively. Note that
$X,Y$ and $Z$ satisfy stochastic differential equations

$\left\\{\begin{array}[]{l}dX=\gamma_{1}X[(1-X)(X-\alpha_{1})+\beta_{12}(Y-y_{0})+\beta_{13}(Z-z_{0})+\frac{\sigma_{1}^{2}}{2}]dt+\sigma_{1}XdW_{1};\\\
dY=\gamma_{2}Y[y_{0}-Y+\beta_{21}X+\beta_{23}(Z-z_{0})+\frac{\sigma_{2}^{2}}{2}]dt+\sigma_{2}YdW_{2};\\\
dZ=\gamma_{3}Z[(z_{0}-Z)(Z-\alpha_{2})+\beta_{31}X+\beta_{32}(Y-y_{0})+\frac{\sigma_{3}^{2}}{2}]dt+\sigma_{3}ZdW_{3};\\\
X(0)=x(0),~{}~{}~{}~{}~{}Y(0)=y(0),~{}~{}~{}~{}~{}Z(0)=z(0).\end{array}\right.$
(11)

Derivation of (11) can be found in [8] where it is stated as Theorem 4.2.1.
Here and below, capital letters denote random variables and processes.

In what follows, all the coefficients in (7), (9), (10) and (11) are assumed
to be positive apart from $\beta_{12}\leq 0$ and $\beta_{32}\leq 0$.

######  Lemma 4.1.

Stochastic differential equations (10) (and hence (11)) have a unique strong
continuous solution on the time horizon $[0,\infty)$.

The proof is presented in the Appendix. It implies that, under positive
initial conditions, the ordinary differential equations (7) (and hence (6))
have a unique solution such that $x(t),y(t),z(t)>0$ for all $t\geq 0$: see
Remark 7.2.

Figure 7:  Stochastic dynamics of the ethnogenesis: examples of the solution
to stochastic differential equations (11). Parameters values:
$\alpha_{1}=0.03$, $\alpha_{2}=0.11$, $y_{0}=0.075$, $z_{0}=0.22$,
$\beta_{12}=-6$, $\beta_{13}=0.6$, $\beta_{21}=0.2$, $\beta_{23}=0.1$,
$\beta_{31}=0.5$, $\beta_{32}=0$, $\gamma_{1}=1$, $\gamma_{2}=0.7$,
$\gamma_{3}=0.2$, Initial conditions: $X(0)=0.07$, $Y(0)=0.053$, $Z(0)=0.05$.
Panel A: $\sigma_{1}=\sigma_{2}=\sigma_{3}=0.05$. Panel B:
$\sigma_{1}=\sigma_{2}=\sigma_{3}=0.1$.

Examples of the stochastic dynamics exhibited in the system (11) are shown in
Fig. 7. The values of model parameters used for this illustration are
identical to those in Fig. 4. The amplitude of noise in Panel B is twice
higher than in Panel A. Solid lines in the both panels show the stochastic
dynamics for three subgroups composing the population. For comparison, we also
provide the deterministic curves which are represented by the dashed lines
(note, that they are identical to the lines shown in Fig. 4, panel B). It is
evident that the stochastic dynamics is significantly different from the
deterministic one. One can see from Panel A that the accumulation of noise
results in the much higher bust in the level of passionaries: the amplitude of
the bust in the stochastic case is about 0.7 against 0.5 in the deterministic
one. Also in Panel A we see that the noise causes significant change in the
level of subpassionaries: at t=100 this level in stochastic case is over 0.3
while it should be about 0.13 in the deterministic case. The impact of noise
is even more evident from Panel B where the amplitude of noise is twice higher
than in Panel A. We can see that the increase in the level of noise not only
increases the discrepancy with the deterministic case (the amplitude of bust
in Panel B is over 0.8) but also results in the occurrence of a new bust.
While the first bust in Panel B was initiated manually, the second bust
appears due to the stochastic effects in the system. This observation leads us
to the conclusion that the ethnogenesis can be initiated by the noise in the
environment surrounding the population.

##  5 Interaction of Ethnogenetic Processes

We conclude our study with modelling the interaction between two ethnoses,
following the same ethnogenetic processes, which however are shifted over
time. The (random) sizes of subpopulations of passionaries, harmonious people
and subpassionaries for the first ethnos are denoted as $X_{1}$, $Y_{1}$ and
$Z_{1}$, while for the second ethnos as $X_{2}$, $Y_{2}$ and $Z_{2}$. We
assume that, being isolated, the ethnoses are identical, described by the
stochastic differential equations like (11), but influenced by six mutually
independent Brownian motions $W_{1},W_{2},\ldots,W_{6}$. We introduce the time
lag between two ethnogenetic processes, such that the first one starts at
$T=0$, the second ethnos appears $T_{1}$ time units later than the first one,
and communication begins $T_{2}$ time units later, at the time moment
$T_{1}+T_{2}$. For simplicity, we also assume that communication is only among
the passionaries, and they suppress each other.

Therefore, we investigate the following system of six stochastic differential
equations

$\left\\{\begin{array}[]{rcl}dX_{1}&=&\gamma_{1}X_{1}[(1-X_{1})(X_{1}-\alpha_{1})+\beta_{12}(Y_{1}-y_{0})+\beta_{13}(Z_{1}-z_{0})+\frac{\sigma_{1}^{2}}{2}\\\
&&-c_{1}\mathbb{I}\\{t\geq T_{1}+T_{2}\\}X_{2}]dt+\sigma_{1}X_{1}dW_{1};\\\
dY_{1}&=&\gamma_{2}Y_{1}[y_{0}-Y_{1}+\beta_{21}X_{1}+\beta_{23}(Z_{1}-z_{0})+\frac{\sigma_{2}^{2}}{2}]dt+\sigma_{2}Y_{1}dW_{2};\\\
dZ_{1}&=&\gamma_{3}Z_{1}[(z_{0}-Z_{1})(Z_{1}-\alpha_{2})+\beta_{31}X_{1}+\beta_{32}(Y_{1}-y_{0})+\frac{\sigma_{3}^{2}}{2}]dt+\sigma_{3}Z_{1}dW_{3};\\\
X_{1}(0)&=&x(0),~{}~{}~{}~{}~{}Y_{1}(0)=y(0),~{}~{}~{}~{}~{}Z_{1}(0)=z(0);\\\
{}{}\hfil\\\ dX_{2}&=&\mathbb{I}\\{t\geq
T_{1}\\}\Bigl{\\{}\gamma_{1}X_{2}[(1-X_{2})(X_{2}-\alpha_{1})+\beta_{12}(Y_{2}-y_{0})+\beta_{13}(Z_{2}-z_{0})+\frac{\sigma_{1}^{2}}{2}\\\
&&-c_{2}\mathbb{I}\\{t\geq
T_{1}+T_{2}\\}X_{1}]dt+\sigma_{1}X_{2}dW_{4}\Bigr{\\}};\\\
dY_{2}&=&\mathbb{I}\\{t\geq
T_{1}\\}\Bigl{\\{}\gamma_{2}Y_{2}[y_{0}-Y_{2}+\beta_{21}X_{2}+\beta_{23}(Z_{2}-z_{0})+\frac{\sigma_{2}^{2}}{2}]dt+\sigma_{2}Y_{2}dW_{5}\Bigr{\\}};\\\
dZ_{2}&=&\mathbb{I}\\{t\geq
T_{1}\\}\Bigl{\\{}\gamma_{3}Z_{2}[(z_{0}-Z_{2})(Z_{2}-\alpha_{2})+\beta_{31}X_{2}+\beta_{32}(Y_{2}-y_{0})+\frac{\sigma_{3}^{2}}{2}]dt+\sigma_{3}W_{6}\Bigr{\\}};\\\
X_{2}(0)&=&X_{2}(T_{1})=x(0),~{}~{}~{}~{}~{}Y_{2}(0)=Y_{2}(T_{1})=y(0),~{}~{}~{}~{}~{}Z_{2}(0)=Z_{2}(T_{1})=z(0).\\\
\end{array}\right.$ (12)

The meaning of all the parameters is the same as in the previous models (i.e.
model (11)). Two new parameters ($c_{1},c_{2}>0$) define the strength of
suppressive interactions between passionaries in the two ethnic groups. This
system of stochastic differential equations has a unique strong continuous
solution on the time horizon $[0,\infty)$. The proof of this statement is
similar to the proof of Lemma 4.1.

Figure 8:  Stochastic dynamics of comunicating ethnoses. Parameters values:
$\alpha_{1}=0.03$, $\alpha_{2}=0.11$, $y_{0}=0.075$, $z_{0}=0.22$,
$\beta_{12}=-6$, $\beta_{13}=0.6$, $\beta_{21}=0.2$, $\beta_{23}=0.1$,
$\beta_{31}=0.5$, $\beta_{32}=0$, $\gamma_{1}=1$, $\gamma_{2}=0.7$,
$\gamma_{3}=0.2$, $c_{1}=c_{2}=0.22$. The black vertical lines show the moment
of birth of the younger ethnos and the moment when the communication begins.
Initial conditions: $X(0)=0.07$, $Y(0)=0.053$, $Z(0)=0.05$.

Two examples of dynamics in the interacting ethnoses, described by the system
(12), is given in Fig.8. Only the bust dynamics exhibited by the passionaries
in both ethnoses is shown on this figure, with the solid lines showing
stochastic dynamics, dashed - deterministic
($\sigma_{1}=\sigma_{2}=\sigma_{3}=0$) and dotted - the deterministic dynamics
in the case of non-interacting ethnoses ($c_{1}=c_{2}=0$). Dotted lines have
identical shapes and this indicates that the two ethnogenetic processes, in
the absence of the noise and interaction between the ethnoses, are identical.
While the dashed blue line is almost identical to the dotted blue line, the
dashed red line is considerably lower than the dotted red line, and this
indicates that, in the absence of the noise, the younger ethnos (dashed red
line) is suppressed by the older ethnos (dashed blue line). Finally we note
that the solid blue line in Panel A is higher than the dashed blue line, while
the solid red line is lower than the dashed red line. This observation
illustrated the impact of the noise to the dynamics of the interacting
ethnoses, which in this particular case results in the amplification of the
suppression of the younger ethnos by the older one. Looking at the shapes of
the solid and dashed lines on Panel B we come to the conclusion that the noise
can also result in the suppression of the older ethnos by the younger one.
Comparing dynamics presented in Panels A and B we note that the dynamics
exhibited by two interacting ethnoses is greatly affected by the noise,
although, as numerous simulations confirm, the scenario from Panel A is more
likely to take place.

##  6 Discussion

In this work we have presented the mathematical model of ethnogenesis which we
have developed on the basis of the paradigm of ”passionary tension” introduced
by Gumilev [4]. According to Gumilev, passionary tension can occur in certain
polities as a result of formation and growth of a subgroup of positively
motivated people, whom Gumilev called ”passionaries”. The idea that the growth
and evolution of a polity (which can cause the formation of civilisation) is
based on its internal structure, and particularly, on the formation of a
certain subgroup of people who push the polity forward, was introduced by
arabic historian Ibn Khaldun in the 15th century [6]. One of the main points
made by Gumilev is that the measure of passionary tension in the polity is
given by its size, i.e. the size of population or territory. As the
quantitative data on the territorial expansion and collapse of past
civilisations are known much better than on their population sizes, it makes
sense to use the size of area taken by a polity as a measure of the passionary
tension in this polity.

The model we have presented here is based on the consideration of the internal
structure of the polity with the dynamics of this structure described by
ordinary differential equations. The main point about the polity’s internal
structure is that there is a subgroup of people, namely, passionaries, and the
size of this subgroup gives a measure of the passionary tension in the polity,
which in turn can be considered as the measure of the size of territory
occupied by the polity. This approach allows to consider the interaction of
the given polity with its neighbours indirectly: the polity’s geopolitical
success is proportional to the number of passionaries in it.

The main feature of the model we have presented here is that it produces the
excitable dynamics in the structure of the evolving polity. That is, when the
polity is in equilibrium (in homoeostatic state) there are no passionaries in
it. However, if there appear a small number of passionaries, this number grows
up to considerable level and then declines back to zero. Thus, formation of
busts, describing the raise and fall of civilisations, is considered here as
an excitation process. Using different versions of the model we have performed
the following studies:

  * •

In the two-variable model given by (6) we explored the types of interactions
between the passionaries and the rest of the population resulting in the
excitable kinetics. Two types of such interactions are illustrated in Figures
2 and 3.

  * •

In the three-variable model given by (7) we explored the types of interactions
between three groups allowing to observe the excitable kinetics. Three types
of such interactions are illustrated in Figures 4, 5 and 6.

  * •

In the stochastic model given by (11) we have found that the noise, when
imposed into the model equations, tends to amplify and results in significant
variations in the amplitude of the bust in the system (see Figure 7). One can
conclude that such noise probably adds to the variation of the territorial
size and duration of life of different civilisations.

  * •

In the model of interaction polities given by (12) we studied the interaction
of the polities of different age undergoing ethnogenesis. We found that if
they interact in a way that passionaries from one polity suppress the
passionaries from the other, then the older polity will be more successful, if
the success is measured by the number of passionaries in the polity. However,
this is not necessarily the case when we impose the noise (see Figure 8).

The presented model can be extended in various ways for further studies. One
of such studies can focus on the interaction of polities under a range of
different assumptions about the ways these polities interact. Another obvious
direction for future research is to extend the model in order to fit it to
available observation data.

##  7 Appendix

Proof of Lemma 4.1. Similarly to Fig. 4, panel A, we present the nullclines
$\dot{x}=0$ and $\dot{z}=0$ at $y=0$ in Fig. 9, panel A. Fix a point
$(a,c)\in\mathbb{R}_{+}\times\mathbb{R}_{+}$ such that

$\displaystyle(1-a)(a-\alpha_{1})-\beta_{12}y_{0}+\beta_{13}(z-z_{0})<0$
$\displaystyle\mbox{for all }z\leq c;$
$\displaystyle(z_{0}-c)(c-\alpha_{2})+\beta_{31}x-\beta_{32}y_{0}<0$
$\displaystyle\mbox{for all }x\leq a.$

The half-open rectangle $(0,a]\times(0,c]$ is shown with the green lines.
Clearly, it is always possible to increase simultaneously $a$ and $c$, so we
assume that $a,c>1$, $0<x(0)<a$ and $0<z(0)<c$. Roughly speaking, the point
$(a,c)$ is outside the ‘internal part’ of the both parabolas.

The similar picture in the variables (8) is given in Fig. 9, panel B: the
images of the parabolas, shown with the blue lines, represent the nullclines
$\dot{v}_{1}=0$ and $\dot{v}_{3}=0$ of equation (9) in the limiting case when
$v_{2}\to-\infty$.

Figure 9: Panel A: nullclines $\dot{x}=0$ and $\dot{z}=0$ at $y=0$ are shown
in blue. Panel B: nullclines $\dot{v}_{1}=0$ and $\dot{v}_{3}=0$ at
$v_{2}\to-\infty$ are shown in blue. The dotted blue line corresponds to the
part of the parabola $\dot{z}=0$ with negative values of $x$. Parameters
values: $\alpha_{1}=0.03$, $\alpha_{2}=0.11$, $y_{0}=0.075$, $z_{0}=0.22$,
$\beta_{12}=-6$, $\beta_{13}=0.6$, $\beta_{21}=0.2$, $\beta_{23}=0.1$,
$\beta_{31}=0.5$, $\beta_{32}=0$, $\gamma_{1}=1$, $\gamma_{2}=0.7$,
$\gamma_{3}=0.2$.

In the space $\mathbb{R}^{3}$, consider the infinite closed prism $\Pi$
defined by

$-\infty<v_{1}\leq\ln a,~{}~{}~{}~{}~{}-\infty<v_{2}\leq\ln
b,~{}~{}~{}~{}~{}-\infty<v_{3}\leq\ln c,$

where

$b>y_{0}+\beta_{21}a+\beta_{23}(c-z_{0}).$ (13)

Without loss of generality, we assume that $y(0)<b$ and $b>1$.

######  Remark 7.1.

For such a prism, we have the following.

  * •

If $x=a$, then the square bracket in the first equation (7) is negative for
all $z\leq c$ and all $y>0$. (Recall that $\beta_{12}\leq 0$.) Therefore, if
$v_{1}=\ln a$, then the square bracket in the first equation (9) and (10) is
negative for all $v_{3}\leq\ln c$ and all $v_{2}>-\infty$.

  * •

Similarly, if $z=c$ then the square bracket in the third equation (7) is
negative for all $x\leq a$ and all $y>0$. (Recall that $\beta_{32}\leq 0$.)
Therefore, if $v_{3}=\ln c$, then the square bracket in the third equation (9)
and (10) is negative for all $v_{1}\leq\ln a$ and all $v_{2}>-\infty$.

  * •

If $y=b$, then the square bracket in the second equation (7) is negative for
all $x\leq a$ and $z\leq c$. Therefore, if $v_{2}=\ln b$, then the square
bracket in the second equation (9) and (10) is negative for all $v_{1}\leq\ln
a$ and all $v_{3}\leq\ln c$.

Now we modify the equations (10) outside $\Pi$: if
$(v_{1},v_{2},v_{3})=m(\hat{v}_{1},\hat{v}_{2},\hat{v}_{3})$ for some $m>1$
with $(\hat{v}_{1},\hat{v}_{2},\hat{v}_{3})\in\partial\Pi$, then we put

$\left\\{\begin{array}[]{l}f_{1}(v_{1},v_{2},v_{3}):=\gamma_{1}[(1-e^{\hat{v}_{1}})(e^{\hat{v}_{1}}-\alpha_{1})+\beta_{12}(e^{\hat{v}_{2}}-y_{0})+\beta_{13}(e^{\hat{v}_{3}}-z_{0})];\\\
f_{2}(v_{1},v_{2},v_{3}):=\gamma_{2}[y_{0}-e^{\hat{v}_{2}}+\beta_{21}e^{\hat{v}_{1}}+\beta_{23}(e^{\hat{v}_{3}}-z_{0})]dt;\\\
f_{3}(v_{1},v_{2},v_{3}):=\gamma_{3}[(z_{0}-e^{\hat{v}_{3}})(e^{\hat{v}_{3}}-\alpha_{2})+\beta_{31}e^{\hat{v}_{1}}+\beta_{32}(e^{\hat{v}_{2}}-y_{0})]dt\end{array}\right.$

and introduce stochastic differential equations (further, SDEs)

$\left\\{\begin{array}[]{l}dV_{1}=f_{1}(V_{1},V_{2},V_{3})dt+\sigma_{1}dW_{1};\\\
dV_{2}=f_{2}(V_{1},V_{2},V_{3})dt+\sigma_{2}dW_{2};\\\
dV_{3}=f_{3}(V_{1},V_{2},V_{3})dt+\sigma_{3}dW_{3},\\\
V_{1}(0)=v_{1}(0),~{}~{}~{}~{}~{}V_{2}(0)=v_{2}(0),~{}~{}~{}~{}~{}V_{3}(0)=v_{3}(0).\end{array}\right.$
(14)

They satisfy all the conditions which guarantee the existence of the unique
continuous strong solution [1, Remark 14.21] or [8, Theorem 5.2.1]: all the
functions $f_{1},f_{2}$ and $f_{3}$ are bounded and Lipschitz in
$\mathbb{R}^{3}$.

######  Remark 7.2.

If $\sigma_{1}=\sigma_{2}=\sigma_{3}=0$, we have just the system of ordinary
differential equations which is uniquely solvable on the time horizon
$[0,\infty)$ [9, Corollary 2.6.]. This solution $(v_{1}(t),v_{2}(t),v_{3}(t))$
can never leave the prism $\Pi$ because on the bounday $\ln a$ for the
component $v_{1}$ ($\ln b$ for $v_{2}$ and $\ln c$ for $v_{3}$) the derivative
$\dot{v}_{1}$ is negative ($\dot{v}_{2}<0$ and $\dot{v}_{3}<0$
correspondingly): see Remark 7.1.

Finally, within the prism $\Pi$, the vector $(v_{1}(t),v_{2}(t),v_{3}(t))$
satisfies differential equations (9) and the functions $x(t):=e^{v_{1}(t)}$,
$y(t):=e^{v_{2}(t)}$ and $z:=e^{v_{3}(t)}$ are well defined on the infinite
horizon $[0,\infty)$, satisfy equations (7) and are strictly positive. As was
noted below (7), these equations cannot have other solutions.

In the stochastic version with $\sigma_{1},\sigma_{2},\sigma_{3}>0$, the
solution to SDE (14) can exit any one prism on a finite time interval. We need
to define a sequence of increasing prisms $\\{\Pi_{i}\\}_{i=0}^{\infty}$
coming from a carefully selected sequence
$\\{(a_{i},b_{i},c_{i})\\}_{i=0}^{\infty}$. Namely, we require that, for a
preliminarily fixed $k>0$, the following condition is satisfied.

######  Condition 7.1.

For each $i>0$ for all

$\displaystyle a_{i-1}\leq x\leq a_{i},~{}~{}~{}~{}~{}0<y\leq
b_{i},~{}~{}~{}~{}~{}0<z\leq c_{i}$ $\displaystyle(\mbox{or}$ $\displaystyle
0<x\leq a_{i},~{}~{}~{}~{}~{}b_{i-1}\leq y\leq b_{i},~{}~{}~{}~{}~{}0<z\leq
c_{i}$ or $\displaystyle 0<x\leq a_{i},~{}~{}~{}~{}~{}0<y\leq
b_{i},~{}~{}~{}~{}~{}c_{i-1}\leq z\leq c_{i})$

the square bracket in the first (second, third) equation (7) is negative and
$a_{i}\geq a_{i-1}e^{k}$ ($b_{i}\geq b_{i-1}e^{k}$, $c_{i}\geq c_{i-1}e^{k}$
correspondingly). As the result, for all $(v_{1},v_{2},v_{3})\in\Pi_{i}$ with
$v_{1}\in[\ln a_{i-1},\ln a_{i}]$ (with $v_{2}\in[\ln b_{i-1},\ln b_{i}]$,
$v_{3}\in[\ln c_{i-1},\ln c_{i}]$) the square bracket in the first
(correspondingly, second, third) equation (9) and (10) is negative.

Additionally, $\ln a_{i}\geq\ln a_{i-1}+k$ ($\ln b_{i}\geq\ln b_{i-1}+k$, $\ln
c_{i}\geq\ln c_{i-1}+k$) and $x(0)<a_{0},y(0)<b_{0},z(0)<c_{0}$.

Clearly, under this condition,

$\inf\\{|\vec{u}_{i}-\vec{u}_{i+1}|:~{}~{}\vec{u}_{i}\in\partial\Pi_{i},~{}\vec{u}_{i+1}\in\partial\Pi_{i+1}\\}\geq
k>0.$

Let us explain why Condition 7.1 can be satisfied for an arbitrarily fixed
$k>0$.

Along with the parabolas as in Fig. 9, panel A, we introduce the expanded
graphs (shown in Fig. 10, panel A with the dashed blue lines) of the functions

$\displaystyle z$ $\displaystyle=$
$\displaystyle\frac{-1}{\beta_{13}e^{k}}\left[(1-x)(x-\alpha_{1})-\beta_{12}y_{0}-\beta_{13}z_{0}\right];$
$\displaystyle x$ $\displaystyle=$
$\displaystyle\frac{-1}{\beta_{31}e^{k}}\left[(z_{0}-z)(z-\alpha_{2})-\beta_{32}y_{0}\right].$

The right-hand parts become bigger than $z$ and $x$ for big enough $x$ and $z$
correspondingly, and one can choose $a_{0}=c_{0}$ such that

$\displaystyle a_{0}=c_{0}$ $\displaystyle<$
$\displaystyle\frac{-1}{\beta_{13}e^{k}}\left[(1-x)(x-\alpha_{1})-\beta_{12}y_{0}-\beta_{13}z_{0}\right];$
$\displaystyle a_{0}=c_{0}$ $\displaystyle<$
$\displaystyle\frac{-1}{\beta_{31}e^{k}}\left[(z_{0}-z)(z-\alpha_{2})-\beta_{32}y_{0}\right]$
and $\displaystyle a_{0}>x(0),~{}~{}~{}~{}~{}c_{0}>z(0).$

After that, the whole red square in Fig. 10, panel A, with $a_{1}=a_{0}e^{k}$
and $c_{1}=c_{0}e^{k}=a_{1}$, is within the area where the square brackets in
the first and third equations (7) are negative. The image of Fig. 10, panel A
on the plain $(v_{1},v_{3})$ is given in Fig. 10, panel B. It remains to take

$b_{0}>\max\\{y_{0}+\beta_{21}a_{1}+\beta_{23}(c_{1}-z_{0}),~{}~{}~{}y(0)\\}.$

In general, for all $i\geq 1$, we put

$\displaystyle a_{i}=a_{i-1}e^{k},~{}~{}~{}~{}~{}c_{i}=c_{i-1}e^{k}$ and
$\displaystyle
b_{i}=\max\\{y_{0}+\beta_{21}a_{i+1}+\beta_{23}(c_{i+1}-z_{0}),~{}~{}b_{i-1}e^{k}\\}.$

The obtained sequence $\\{(a_{i},b_{i},c_{i})\\}_{i=0}^{\infty}$ satisfies
Condition 7.1.

Figure 10: Construction of the first prisms $\Pi_{0}$ and $\Pi_{1}$. Parameter
values are as in Fig. 9; $k=0.75$.

The SDE (14), for the prism $\Pi_{i}$, in its vector form means that

$\vec{V}^{i}(t,\omega)=\vec{v}(0)+\int_{0}^{t}\vec{f}^{i}(\vec{V}^{i}(s,\omega))ds+\Xi\vec{W}(t,\omega)~{}~{}~{}~{}~{}a.s.$
(15)

for all $t\geq 0$. Here the functions
$\vec{f}^{i}=(f^{i}_{1},f^{i}_{2},f^{i}_{3})$ are constructed, as described
above, for the prism $\Pi_{i}$. The vector notations are conventional,
$\Xi=\left[\begin{array}[]{ccc}\sigma_{1}&0&0\\\ 0&\sigma_{2}&0\\\
0&0&\sigma_{3}\end{array}\right]$. Similarly, we write down equations (10) as

$\vec{V}(t,\omega)=\vec{v}(0)+\int_{0}^{t}\vec{f}(\vec{V}(s,\omega))ds+\Xi\vec{W}(t,\omega)~{}~{}~{}~{}~{}a.s.$
(16)

for all $t\geq 0$. Note that, if a random process $\vec{Z}(t,\omega)$
satisfies equation (15), then $\vec{V}^{i}(t,\omega)=\vec{Z}(t,\omega)$ for
all $t\geq 0$ almost surely. (The processes $\vec{V}^{i}$ and $\vec{Z}$ are
‘indistinguishable’, often called ’versions’ or ’modifications’ [8].)

Suppose $\tau>0$ is arbitrarily fixed, construct the process $\vec{V}$ on
$[0,\tau]$, which satisfies equation (16), and prove that it is unique. The
idea is as follows.

  * •

The process $\vec{V}$ will be a combination of the processes $\vec{V}^{i}$.

  * •

Between the prisms $\Pi_{i-1}$ and $\Pi_{i}$, the process $\vec{V}^{i}$ is
pushed back to the prism $\Pi_{i-1}$, and the chance for it to leave the prism
$\Pi_{i-1}$ is smaller than a constant $\varepsilon<1$.

  * •

Hence, almost surely, there is $j\geq 0$ such that the process $\vec{V}^{j}$
lives in the prism $\Pi_{j}$, where it is unique and coincides with $\vec{V}$.

After that, one can extend the unique solution to (16) to the infinite horizon
$[0,\infty)$.

We introduce ‘debutes’

$D_{i}(\omega):=\inf\\{0<t\leq\tau:~{}~{}\vec{V}^{i}(t,\omega)\in\Pi_{i+1}\setminus\Pi_{i}\\},~{}~{}~{}i=0,1,\ldots$

and the (measurable) sets

$\Omega_{i}:=\\{\omega\in\Omega:~{}~{}D_{i}(\omega)\geq\tau\\},~{}~{}i=0,1,\ldots.$

Any debute is a Markov moment because all the processes $\vec{V}^{i}$ are
continuous a.s. [8, Example 7.2.2.]. On each set $\Omega_{i}$, for all
$s\leq\tau$, $\vec{V}^{i}(s,\omega)\in\Pi_{i}$ meaning that
$\vec{f}^{i}(\vec{V}^{i}(s,\omega))=\vec{f}(\vec{V}^{i}(s,\omega))$ and, by
(15), for all $t\leq\tau$, for $P$-almost all $\omega\in\Omega_{i}$,

$\vec{V}^{i}(t,\omega)=\vec{v}(0)+\int_{0}^{t}\vec{f}(\vec{V}^{i}(s,\omega))ds+\Xi\vec{W}(t,\omega).$
(17)

For $0\leq s\leq\tau$, we put

$\displaystyle\vec{V}(s,\omega)$ $\displaystyle:=$
$\displaystyle\vec{V}^{0}(s,\omega)~{}\mbox{ on
}\Omega_{0}~{}~{}~{}~{}~{}\mbox{ and }$ $\displaystyle\vec{V}(s,\omega)$
$\displaystyle:=$ $\displaystyle\vec{V}^{i}(s,\omega)~{}\mbox{ on
}\overline{\cup_{j=0}^{i-1}\Omega_{j}}\cap\Omega_{i}~{}~{}~{}~{}~{}\mbox{ for
}i=1,2,\ldots.$

As the result, the continuos process $\vec{V}$ is built on
$\cup_{i=0}^{\infty}\Omega_{i}$.

Let us show that

$P\left(\bigcup_{i=0}^{\infty}\Omega_{i}\right)=1.$ (18)

This will be done in three steps.

In what follows, we use expression ‘on $\Omega^{\prime}\subset\Omega$
statement$(\omega)$ holds a.s.’ to say that
$P(\\{\omega\in\Omega^{\prime}:~{}\mbox{statement}(\omega)~{}\mbox{is
false}\\})=0$.

Step 1. For $i\geq 1$, let us also introduce debutes

$D_{i,i-1}(\omega):=\inf\\{0<t\leq\tau:~{}\vec{V}^{i}(t,\omega)\in\Pi_{i}\setminus\Pi_{i-1}\\}$

and the corresponding subsets

$\Omega_{i,i-1}:=\\{\omega\in\Omega:~{}~{}D_{i,i-1}(\omega)\geq\tau\\},$

and show that

$P(\Omega_{i-1}\triangle\Omega_{i,i-1})=0.$ (19)

Consider the modified $\vec{V}^{i-1}$ process

$\vec{V}^{i-1^{\prime}}(t,\omega):=\left\\{\begin{array}[]{ll}\vec{V}^{i}(t,\omega),&\mbox{
if }\omega\in\Omega_{i,i-1};\\\ \vec{V}^{i-1}(t,\omega),&\mbox{ if
}\omega\in\overline{\Omega_{i,i-1}}\end{array}\right.$

and show that

$\vec{V}^{i-1^{\prime}}(t,\omega)=\vec{V}^{i-1}(t,\omega)~{}~{}~{}~{}~{}a.s.$
(20)

for all $t\in[0,\tau]$.

On the set $\overline{\Omega_{i,i-1}}$, we have, for all $t\in[0,\tau]$,

$\vec{V}^{i-1^{\prime}}(t,\omega)=\vec{v}(0)+\int_{0}^{t}\vec{f}^{i-1}(\vec{V}^{i-1^{\prime}}(s,\omega))ds+\Xi\vec{W}(t,\omega)~{}~{}~{}~{}~{}a.s.$
(21)

because here $\vec{V}^{i-1^{\prime}}=\vec{V}^{i-1}$. On the set
$\Omega_{i,i-1}$, $\vec{V}^{i-1^{\prime}}=\vec{V}^{i}\in\Pi_{i-1}$ for all
$t\in[0,\tau]$, and, for the vectors $\vec{v}$ from $\Pi_{i-1}\subset\Pi_{i}$
we have $\vec{f}^{i-1}(\vec{v})=\vec{f}^{i}(\vec{v})=\vec{f}(\vec{v})$ meaning
that for all $t\in[0,\tau]$ again

$\vec{V}^{i-1^{\prime}}(t,\omega)=\vec{v}(0)+\int_{0}^{t}\vec{f}^{i-1}(\vec{V}^{i-1^{\prime}}(s,\omega))ds+\Xi\vec{W}(t,\omega)~{}~{}~{}~{}~{}a.s.$
(22)

From (21) and (22) we deduce that $\vec{V}^{i-1^{\prime}}$ is a solution to
the SDE (15) at $i-1$, and assertion (20) follows. Therefore, on
$\Omega_{i,i-1}$, $\vec{V}^{i-1}=\vec{V}^{i}\in\Pi_{i-1}$ for all
$t\in[0,\tau]$ almost surely, and $P(\Omega_{i-1}\setminus\Omega_{i,i-1})=0$.

In the similar way, we consider the modified $\vec{V}^{i}$ process

$\vec{V}^{i^{\prime}}(t,\omega):=\left\\{\begin{array}[]{ll}\vec{V}^{i-1}(t,\omega),&\mbox{
if }\omega\in\Omega_{i-1};\\\ \vec{V}^{i}(t,\omega),&\mbox{ if
}\omega\in\overline{\Omega_{i-1}},\end{array}\right.$

which, for all $t\in[0,\tau]$, satisfies the SDE (15) at $i$. Therefore, on
$\Omega_{i-1}$, $\vec{V}^{i-1}=\vec{V}^{i}\in\Pi_{i-1}$ for all $t\in[0,\tau]$
almost surely, and $P(\Omega_{i,i-1}\setminus\Omega_{i-1})=0$. Equality (19)
is proved.

Step 2. Suppose $k$ is big enough. In fact, the choice of $k$ depends on
$\tau$ and $\max\\{\sigma_{1},\sigma_{2},\sigma_{3}\\}$ only: $k$ must only
satisfy inequality

$\frac{4}{\sqrt{2\pi}}\int_{d_{i}}^{\infty}e^{-\frac{y^{2}}{2}}dy<\frac{1}{3},~{}~{}~{}~{}~{}i=1,2,3,$
(23)

where $d_{i}=\frac{k}{2\sigma_{i}\sqrt{\tau}}$.

We are going to show that, for each $i\geq 1$, which is fixed below,

$P(\overline{\Omega_{i}})\leq\varepsilon P(\overline{\Omega_{i-1}}),$ (24)

where, $\varepsilon<1$ is some $i$-independent constant.

Clearly, $\overline{\Omega_{i}}\subset\overline{\Omega_{i,i-1}}$ and, for
$\omega\in\overline{\Omega_{i}}$, $0<D_{i,i-1}(\omega)<D_{i}(\omega)<\tau$
because $\Pi_{i-1}\subset\Pi_{i}$ and the process $\vec{V}^{i}$ is continuous.
We will estimate $P(\overline{\Omega_{i}}|\overline{\Omega_{i,i-1}})$ assuming
that $P(\overline{\Omega_{i,i-1}})>0$. (Otherwise, inequality (24) is
trivial.)

Here and below, usually, all the statements hold $P$-a.s., and all the
introduced random variables are defined for $P$-almost all $\omega\in\Omega$,
without special remarks.

For $\omega\in\overline{\Omega_{i,i-1}}$, the Markov moment
$D_{i,i-1}(\omega)$ is smaller than $\tau$, and we denote
$\vec{U}(\omega):=\vec{V}^{i}(D_{i,i-1}(\omega),\omega)\in\partial\Pi_{i-1}$.
According to the strong Markov property of $\vec{V}^{i}(t,\omega)$ [8, Theorem
7.2.4.], the future behaviour of $\vec{V}^{i}(t,\omega)$ on
$[D_{i,i-1}(\omega),\tau]$ depends only on $\vec{U}(\omega)$. Let us estimate
$P(\overline{\Omega_{i}}|\vec{U}(\omega))$. The set $\overline{\Omega_{i}}$ is
split in three disjoint subsets depending on which component first reaches the
boundary $\partial\Pi_{i}$:

$\displaystyle E_{1}$ $\displaystyle=$
$\displaystyle\\{\omega:~{}V^{i}_{1}(D_{i}(\omega),\omega)=\ln a_{i}\\},$
$\displaystyle E_{2}$ $\displaystyle=$
$\displaystyle\\{\omega:~{}V^{i}_{2}(D_{i}(\omega),\omega)=\ln
b_{i}\\}\setminus E_{1},$ $\displaystyle\mbox{and}~{}~{}E_{3}$
$\displaystyle=$
$\displaystyle\\{\omega:~{}V^{i}_{3}(D_{i}(\omega),\omega)=\ln
c_{i}\\}\setminus(E_{1}\cup E_{2}).$

We shall prove that, for some $i$-independent constant $\delta<\frac{1}{3}$,

$P(E_{j}|U(\omega))\leq\delta,~{}~{}~{}~{}~{}j=1,2,3.$ (25)

Suppose $j=1$: the reasoning for $j=2$ and $j=3$ is similar.

Let $\vec{u}\in\partial\Pi_{i-1}$ be fixed and consider the process
$\vec{V}^{i}(t,\omega)$ with $\omega\in\overline{\Omega_{i,i-1}}$, $t\geq
D_{i,i-1}(\omega)$ as starting from $\vec{V}^{i}(D_{i,i-1}(\omega))=\vec{u}$:

$\vec{V}^{i}(t,\omega)=\vec{u}+\int_{D_{i,i-1}(\omega)}^{t}\vec{f}^{i}(\vec{V}^{i}(s,\omega))ds+\Xi[\vec{W}(t,\omega)-\vec{W}(D_{i,i-1}(\omega),\omega)],~{}t\in[D_{i,i-1}(\omega),\tau].$
(26)

Below, we assume that $\omega\in E_{1}$ in order to estimate
$P(E_{1}|\vec{u})$, so that $D_{i}(\omega)<\tau$ is well defined. The set

$H_{i}(\omega):=\\{t\in[D_{i,i-1}(\omega),D_{i}(\omega)]:~{}V^{i}_{1}(t,\omega)\leq\ln
a_{i-1}\\}$

is not empty: $V^{i}_{1}(D_{i,i-1}(\omega))=u_{1}\leq\ln a_{i-1}$. We put

$T_{i}(\omega):=\sup\\{t:~{}t\in H_{i}(\omega)\\}.$

Although $T_{i}$ is not a Markov moment, it is a measurable random variable:
for each $x\in\mathbb{R}_{+}$

$\displaystyle\\{\omega:~{}T_{i}(\omega)>x\\}$ $\displaystyle=$
$\displaystyle\\{\omega:~{}x<D_{i}(\omega)\mbox{ and }\exists
s\in(x,D_{i}(\omega))\cap\mathbb{Q}:~{}V^{i}_{1}(s,\omega)\leq a_{i-1}\\}$
$\displaystyle=$
$\displaystyle\\{\omega:~{}D_{i}(\omega)>x\\}\cap\bigcup_{s\in\mathbb{Q},~{}s>x}\\{\omega:~{}D_{i}(\omega)>s\mbox{
and }V^{i}_{1}(s,\omega)\leq\ln a_{i-1}\\},$

where $\mathbb{Q}$ is the set of rational numbers. Note that
$V^{i}_{1}(T_{i}(\omega),\omega)=\ln a_{i-1}$ and remember that
$T_{i}(\omega)\in[D_{i,i-1}(\omega),D_{i}(\omega)]$ and
$0<D_{i,i-1}(\omega)<D_{i}(\omega)<\tau$. On the time interval
$[T_{i}(\omega),D_{i}(\omega)]$, the process $\vec{V}^{i}(t,\omega)$ is still
in $\Pi_{i}$ and $V^{i}_{1}(t,\omega)\in[\ln a_{i-1},\ln a_{i}]$: see Fig. 11.

Figure 11: Small fragment of the motion of $\vec{V}^{i}$: phase trajectories
on panel A and time trajectories on Panel B. The black blobs correspond to the
time moments $D_{i,i-1}$, $T$ and $D_{i}$. The red arrows at the blobs on
panel A show the direction of the movement of the phase trajectories at the
time moments $D_{i,i-1}$, $T_{i}$ and $D_{i}$. Parameters values:
$\alpha_{1}=0.03$, $\alpha_{2}=0.11$, $y_{0}=0.075$, $z_{0}=0.22$,
$\beta_{12}=-6$, $\beta_{13}=0.6$, $\beta_{21}=0.2$, $\beta_{23}=0.1$,
$\beta_{31}=0.5$, $\beta_{32}=0$, $\gamma_{1}=1$, $\gamma_{2}=0.7$,
$\gamma_{3}=0.2$, $\sigma_{1}=\sigma_{2}=\sigma_{3}=1$.

According to (26) and taking into account Condition 7.1, we have

$\displaystyle V^{i}_{1}(D_{i}(\omega),\omega)$ $\displaystyle=$
$\displaystyle
u_{1}+\int_{D_{i,i-1}(\omega)}^{T_{i}(\omega)}f^{i}_{1}(\vec{V}^{i}(s,\omega))ds+\sigma_{1}[W_{1}(T_{i}(\omega),\omega)-W_{1}(D_{i,i-1}(\omega),\omega)]$
$\displaystyle+\int_{T_{i}(\omega)}^{D_{i}(\omega)}f^{i}_{1}(\vec{V}^{i}(s,\omega))ds+\sigma_{1}[W_{1}(D_{i}(\omega),\omega)-W_{1}(T_{i}(\omega),\omega)]$
$\displaystyle=$ $\displaystyle
V^{i}_{1}(T_{i}(\omega),\omega)+\int_{T_{i}(\omega)}^{D_{i}(\omega)}f^{i}_{1}(\vec{V}^{i}(s,\omega))ds+\sigma_{1}[W_{1}(D_{i}(\omega),\omega)-W_{1}(T_{i}(\omega),\omega)]$
$\displaystyle\leq$ $\displaystyle
V^{i}_{1}(T_{i}(\omega),\omega)+\sigma_{1}[W_{1}(D_{i}(\omega),\omega)-W_{1}(T_{i}(\omega),\omega)]:$

the function $f^{i}_{1}(v)$ is negative for $v\in\Pi_{i}$ with $v_{1}\in[\ln
a_{i-1},\ln a_{i}]$. Since $V^{i}_{1}(T_{i}(\omega),\omega)=\ln a_{i-1}$ and
$V^{i}_{1}(D_{i}(\omega),\omega)=\ln a_{i}$, we see that, for $\omega\in
E_{1}$,

$W_{1}(D_{i}(\omega),\omega)-W_{1}(T_{i}(\omega),\omega)\geq\frac{\ln
a_{i}-\ln a_{i-1}}{\sigma_{1}}\geq\frac{k}{\sigma_{1}}.$

That means

$\sup_{0\leq s\leq\tau}W_{1}(s,\omega)-\inf_{0\leq
s\leq\tau}W_{1}(s,\omega)\geq\frac{k}{\sigma_{1}},$

and $P(\tilde{\Omega})$, the probability of the set $\tilde{\Omega}$ of all
the points $\omega$ satisfying this property, is smaller than some
$0<\delta<1/3$ for big enough $k$. To be more specific,

$\tilde{\Omega}\subset\\{\omega:~{}\sup_{0\leq
s\leq\tau}W_{1}(s,\omega)\geq\frac{k}{2\sigma_{1}}\mbox{ or }\inf_{0\leq
s\leq\tau}W_{1}(s,\omega)\leq\frac{-k}{2\sigma_{1}}\\},$

and the probability of the set on the right is well studied: it is smaller
than

$4P(\\{\omega:~{}W_{1}(1,\omega)\geq
d\\})=\frac{4}{\sqrt{2\pi}}\int_{d}^{\infty}e^{-\frac{y^{2}}{2}}dy,$

where $d=\frac{k}{2\sigma_{1}\sqrt{\tau}}$, and approaches zero as
$k\to\infty$. See also Lemma 7.1. Thus, for fixed $\vec{u}$, if $k$ is big
enough (and certainly dependent on
$\max\\{\sigma_{1},\sigma_{2},\sigma_{3}\\}$ and $\tau$ only), then
$P(E_{1}|\vec{u})=\delta<\frac{1}{3}$: inequality (25) is proved.

Applying the similar reasoning to $E_{2}$ and $E_{3}$, we conclude that
$P(\overline{\Omega_{i}}|\vec{u})\leq\varepsilon:=3\delta<1$ and

$P(\overline{\Omega_{i}})=\int_{\overline{\Omega_{i,i-1}}}P(\overline{\Omega_{i}}|\vec{U}(\omega))dP(\omega)\leq\varepsilon
P(\overline{\Omega_{i,i-1}}).$

Inequality (24) now follows from (19).

Step 3. Using (24), we have

$P\left(\bigcap_{i=0}^{\infty}\overline{\Omega_{i}}\right)=\lim_{N\to\infty}P\left(\bigcap_{i=0}^{N}\overline{\Omega_{i}}\right)\leq\lim_{N\to\infty}P(\overline{\Omega_{N}})\leq\lim_{N\to\infty}\varepsilon^{N}P(\overline{\Omega_{0}})=0.$

Therefore,
$P\left(\bigcup_{i=0}^{\infty}\Omega_{i}\right)=1-P\left(\bigcap_{i=0}^{\infty}\overline{\Omega_{i}}\right)=1$,
and equality (18) is proved.

According to (17), the constructed continuous process $\vec{V}(t,\omega)$
satisfies equation (16) (and (10)) for all $t\in[0,\tau]$ a.s., i.e., it is a
strong continuous solution to those SDEs.

Now let us show that, if $\vec{Z}(t,\omega)$ is a strong continuous solution
to the SDE (16), then, for all $t\in[0,\tau]$,
$\vec{Z}(t,\omega)=\vec{V}(t,\omega)$ a.s.

Since the process $\vec{Z}$ is continuous a.s., it is bounded a.s. and hence
(almost surely) there is (a unique) integer $I(\omega)\geq 0$ such that, for
all $t\in[0,\tau]$, $\vec{Z}(t,\omega)\in\Pi_{I(\omega)}$ and, for
$I(\omega)>0$, there is $t\in[0,\tau]$ such that
$\vec{Z}(t,\omega)\in\Pi_{I(\omega)}\setminus\Pi_{I(\omega)-1}$. Recall that
$\bigcup_{i=0}^{\infty}\Pi_{i}=\mathbb{R}^{3}$, so that, the sets
$\tilde{\Omega}_{i}:=\\{\omega:~{}I(\omega)=i\\}$ are such that
$P\left(\Omega\setminus\bigcup_{i=0}^{\infty}\tilde{\Omega}_{i}\right)=0$. On
the set $\tilde{\Omega}_{i}$, the process $\vec{Z}$ satisfies SDE

$\vec{Z}(t,\omega)=\vec{v}(0)+\int_{0}^{t}\vec{f}^{i}(\vec{Z}(s,\omega))ds+\Xi\vec{W}(t,\omega),~{}~{}~{}~{}~{}t\in[0,\tau]$

for $P$-almost all $\omega\in\tilde{\Omega}_{i}$ because, within the prism
$\Pi_{i}$, $\vec{f}(\vec{v})=\vec{f}^{i}(\vec{v})$. We see that
$\vec{Z}(t,\omega)=\vec{V}^{i}(t,\omega)$ for all $t\in[0,\tau]$ a.s. on
$\tilde{\Omega}_{i}$ and, for $P$-almost all $\omega\in\tilde{\Omega}_{i}$,
$\vec{V}^{i}(t,\omega)\in\Pi_{i}$ for all $t\in[0,\tau]$ and, in case $i>0$,
there is $t\in[0,\tau]$ such that
$\vec{V}^{i}(t,\omega)\in\Pi_{i}\setminus\Pi_{i-1}$. The last assertion means
that, for $P$-almost all $\omega\in\tilde{\Omega}_{i}$, $\omega\in\Omega_{i}$
and, in case $i>0$, $\omega\in\overline{\Omega_{i,i-1}}$. (The process
$\vec{V}^{i}$ left $\Pi_{i-1}$ at some moment.) Since
$\overline{\Omega_{i,i-1}}\triangle\overline{\Omega_{i-1}}=\Omega_{i,i-1}\triangle\Omega_{i-1}$,
for $i>0$, according to (19), $\omega\in\overline{\Omega_{i-1}}$ for
$P$-almost all $\omega\in\tilde{\Omega}_{i}$. We conclude that, for $P$-almost
all $\omega\in\tilde{\Omega}_{i}$,
$\omega\in\overline{\Omega_{i-1}}\cap\Omega_{i}$ ($\omega\in\Omega_{0}$ in
case $i=0$) and $\vec{V}^{i}(t,\omega)=\vec{V}(t,\omega)$ for all
$t\in[0,\tau]$ a.s. on $\tilde{\Omega}_{i}$. As a result,
$\vec{Z}(t,\omega)=\vec{V}(t,\omega)$ for all $t\in[0,\tau]$ a.s on
$\tilde{\Omega}_{i}$. Therefore, $\vec{Z}(t,\omega)=\vec{V}(t,\omega)$ for all
$t\in[0,\tau]$ a.s. on $\bigcup_{i=0}^{\infty}\tilde{\Omega}_{i}$, and the
latter set, as explained above, coincides with $\Omega$ up to a set of
$P$-measure zero.

The uniqueness of the strong continuous solution to SDE (16) (hence, (10)) is
proved.

Finally, extension to the infinite horizon $[0,\infty)$ of the solution
$\vec{V}$ to the SDE (16) (and (10)) is trivial. Take an increasing sequence
$\\{\tau_{j}\\}_{j=1}^{\infty}$, $\tau_{j}>0$, with
$\lim_{j\to\infty}\tau_{j}=\infty$, construct the solutions
$\vec{V}^{\tau_{j}}$ to SDE (16) (and (10)) on the intervals $[0,\tau_{j}]$
and put

$\vec{V}(t,\omega):=\vec{V}^{\tau_{j}}(t,\omega)\mathbb{I}\\{t\in[\tau_{j-1},\tau_{j})\\},$

where $\tau_{0}:=0$. Due to the uniqueness of each process
$\vec{V}^{\tau_{j}}$, for each $k=1,2,\ldots$,
$\vec{V}(t,\omega)=\vec{V}^{\tau_{k}}(t,\omega)$ for all $t\in[0,\tau_{k})$
a.s. and hence the process $\vec{V}$ satisfies SDE (16) (and (10)) for all
$t\in[0,\tau_{k})$ a.s. meaning that $\vec{V}$ satisfies SDE (16) (and (10))
for all $t\in[0,\infty)$ a.s.

If there is another process $\vec{Z}$ satisfying this property then, again due
to the uniqueness of $\vec{V}^{\tau_{j}}$,
$\vec{Z}(t,\omega)=\vec{V}^{\tau_{j}}(t,\omega)$ for all
$t\in[\tau_{j-1},\tau_{j})$ a.s., $j=1,2,\ldots$. Hence
$\vec{Z}(t,\omega)=\vec{V}(t,\omega)$ for all $t\in[0,\infty)$ a.s.

The proof is completed.

######  Lemma 7.1.

Let $W(t,\omega)$ be a Brownian motion on a filtered probability space
$(\Omega,{\cal F},({\cal F}_{t})_{t\geq 0},P)$. Then, for $a,\tau>0$

$P(\sup_{0\leq t\leq\tau}W(t,\omega)\geq a\mbox{ or }\inf_{0\leq
t\leq\tau}W(t,\omega)\leq-a)\leq\frac{4}{\sqrt{2\pi}}\int_{a/\sqrt{\tau}}^{\infty}e^{-\frac{y^{2}}{2}}dy.$

Proof. Let $T_{a}(\omega):=\inf\\{t:~{}W(t,\omega)=a\\}$. Then, by the
symmetry,

$P(\sup_{0\leq t\leq\tau}W(t,\omega)\geq a\mbox{ or }\inf_{0\leq
t\leq\tau}W(t,\omega)\leq-a)\leq
P(T_{a}(\omega)\leq\tau)+P(T_{-a}(\omega)\leq\tau)=2P(T_{a}(\omega)\leq\tau).$
(27)

Now, by the reflection principle, $P(T_{a}\leq\tau)=2P(W(\tau,\omega)\geq a)$,
and we continue (27):

$2P(T_{a}(\omega)\leq\tau)=4\int_{a}^{\infty}\frac{1}{\sqrt{2\pi\tau}}e^{-\frac{z^{2}}{2\tau}}dz=\frac{4}{\sqrt{2\pi}}\int_{a/\sqrt{\tau}}^{\infty}e^{-\frac{y^{2}}{2}}dy.$

The proof is completed.

## References

  * [1] R.J.Elliott. Stochastic Calculus and Applications. Springer-Verlag, NY etc, 1982.
  * [2] R.Fitzhugh. Impulses and physiological states in theoretical models of nerve membrane, Biophysical Journal, 1(6), pp. 445-466, 1961.
  * [3] L.Gumilev. Ethnogenesis and the Biosphere. Progress Publishers, Moscow, 1990.
  * [4] L.N.Gumilev and K.P.Ivanov. Ethnic-processes – the approaches in researches. Sotsiologicheskie Issledovaniya, 1, pp.50–57, 1992, (in Russian).
  * [5] A.K. Guz and L.A. Pautova. Global Ethnosociology. Librokom, Moscow, 2019, (in Russian).
  * [6] Ibn Khaldun. The Muqaddimah: An Introduction to History. Pantheon Books, New York 1958.
  * [7] J.D.Murray. Mathematical Biology I. An Introduction. Springer, New York, 2002.
  * [8] B.Oksendal. Stochastic Differential Equations. Springer, Heidelberg, 2013.
  * [9] G.Teschl. Ordinary Differential Equations and Dynamical Systems. AMS, Rhode Island, 2012.
  * [10] P.Turchin. Historical Dynamics. Princeton University Press, New Jercey, 2003.
  * [11] P.J.Wangersky. Lotka-Volterra population models, Annual Review of Ecology and Systematics, 9, pp.189-218, 1978.