Selected Topics of Social Physics: Nonequilibrium Systems ^†

Abstract

This paper is devoted to nonequilibrium systems in the physics approach to social systems. Equilibrium systems have been considered in the recenly published first part of the review. The style of the paper combines the features of a tutorial and a review, which, from one side, makes it simpler to read for nonspecialists aiming at grasping the basics of social physics, and from the other side, describes several rather recent original models containing new ideas that could be of interest to experienced researchers in the field.

1. Introduction

The term “social physics” was introduced by Comte, who gave the definition [1]: “Social physics is that science which occupies itself with social phenomena, considered in the same light as astronomical, physical, chemical, and physiological phenomena, that is to say as being subject to natural and invariable laws, the discovery of which is the special object of its research.” Comte [2] divided sociology into two main fields, or branches: social statics, or the study of the forces that hold society together; and social dynamics, or the study of the causes of social change. For social dynamics, Comte also used the term “social evolutionism”. The term “social physics” was widely used by Quetelet [3].

The whole review consists of two parts. The first part [4] is devoted to static, or equilibrium, phenomena and systems. In the present part of the review, dynamic, or nonequilibrium, phenomena and systems are considered. This implies that the considered social systems have to be characterised by evolution equations describing how the systems develop and what are their evolutionary stable states.

It goes without saying that in a rather short review, it is impossible to cover the whole of social physics, being a highly multidisciplinary science including a variety of many different topics. Applications of mathematical and physical methods to social systems are really extremely ramified. As examples, it is possible to mention the growth and distribution of cities [5,6,7,8,9,10,11,12,13,14,15], dynamics of urban population [16,17,18,19], optimization of city traffic [20,21,22,23,24,25,26], econophysics [27,28,29,30,31,32,33,34,35,36,37,38], the evolution of cooperation, and, more generally, the interplay between cooperation and competition [39,40,41,42,43,44,45,46,47,48], separation of societies into subgroups and communities [49,50,51,52], voter models [50,53,54,55,56], human-like artificial intelligence [57,58,59,60,61,62,63,64,65], criminology [66,67,68,69,70,71,72,73,74], migration [75,76,77], contagion phenomena and the spread of infections [78,79,80], ecological problems [81,82,83], climate change [84,85,86,87], various complex networks [88,89,90,91,92,93,94,95,96,97,98,99,100], and many other topics [101,102,103].

Since it is impossible to embrace the immensity, it is necessary to select some of topics. The choice of topics in the present review is dictated by the main aims that are: first, to give an introduction into the methods of dealing with dynamical systems that are necessary for analysing the behaviour of social systems. Second, to illustrate these methods on quite simple known models. The last, but not the least, in order to make the present review original, it was necessary to consider some novel, fresh models that have not been mentioned in other review articles. It was natural to choose the examples of new models from those fields that are in the frame of the author’s research interests. Thus the main difference of the present article from the previously published surveys is in two points: (i) more detailed exposition of the methods used for treating social systems; (ii) the consideration of new models that have not been studied in the previous reviews.

The exposition starts with a brief survey of dynamical systems and an overview of the main known evolution equations. Then, several recent models are considered allowing for the description of nontrivial dynamical effects. Evolution equations with functional carrying capacity are introduced. Punctuated evolution is described. A new type of evolution equations characterising symbiosis is presented. The peculiarity of constructing the evolution equations for structural self-organised societies is discussed. Some examples of financial markets are treated, emphasising the important role of herding behaviour.

2. Dynamical Social Systems

The evolution of nonequilibrium social systems is described by the time dependence of observable quantities. The time dependence is typically given through differential equations. Sometimes, one considers discrete time, when the evolution equations are given by difference equations. Generally, time is a continuous variable. This review mainly deals with differential equations.

2.1. Dynamical Systems

The evolution of social systems is described by considering the temporal evolution of observable quantities. As examples of observable quantities, it is possible to mention order parameters, the value of total order, and the fractions of populations. The system of differential equations characterising the system evolution, from the mathematical point of view, is a dynamical system which many studies are devoted to (see, e.g., [104,105,106,107,108]). The main definitions necessary for further consideration are briefly surveyed just below.

Let the set of observable quantities of interest be $x_1,x_2,\dots ,x_N$. Each of the observables varies inside the related manifold ${\mathbb{X}}_n$, which is a part of the real-number axis $\mathbb{R}=(-\infty ,+\infty )$. In order not to overload the text by notations, let us just write that $x\in {\mathbb{R}}^N$.

The set,

$$x\equiv \{x_N:n=1,2,\dots ,N\}\in {\mathbb{R}}^N$$

(1)

is the system state. The manifold ${\mathbb{R}}^N$ can be termed phase space. The time arrow is the variation of time t, beginning at an initial moment $t_0$ and increasing to ∞. Without the loss of generality, it is often admissible to set $t_0=0$. So, often time is treated as the variable in the interval $[0,\infty )$.

The evolution law, or the law of motion, is the dependence on time, t, of the variables $x_n$. The time dependent variable,

$$x=x(t,x_0),$$

(2)

with a state $x_0$ given at a moment of time $t_0$, called initial state,

$$x_0=x(t_0,x_0)=\{x_{0n}:n=1,2,\dots ,N\},$$

(3)

is a dynamical state.

The family of evolution equations forms a general dynamical system. Often, one also requires that the states of a dynamical system would satisfy the group composition law,

$$x(t+t^{\prime },x_0)=x(t,x(t^{\prime },x_0)).$$

(4)

The evolution law (2), with the group property (4), consists of the endomorphisms projecting the phase space onto itself. The family of the endomorphisms composes a flow,

$$\{x(t,\dots ):{\mathbb{R}}_{+}\times {\mathbb{R}}^N⟼{\mathbb{R}}^N\}.$$

(5)

The set of all states forms a trajectory,

$$\{x(t,x_0):t\ge t_0\}.$$

(6)

The states at fixed times are called the points of the trajectory.

General dynamical systems in continuous time are typically represented by differential equations in the normal form,

$$\frac{dx}{dt}=f(x,t),$$

(7)

where

$$f(x,t)=\{f_n(x,t):n=1,2,\dots ,N\}.$$

(8)

The differential equations containing higher-order time derivatives can always be reduced to the normal form by means of renotations. The dynamical system dimension is the number of the evolution equations. In the present case, it is N.

If the right-hand side $f(x,t)$ explicitly depends on time, then Equation (7) is termed a non-autonomous equation, while when $f(x,t)=f\left(x\right)$ does not depend on time, this is an autonomous equation. The intermediate case is when at finite times the evolution equation is non-autonomous, but, as time increases, the equation becomes autonomous, so that

$$\underset{t\to \infty }{lim}f(x,t)=f\left(x\right).$$

Then the equation is called “transiently non-autonomous”. This case corresponds to a nonequilibrium system subject to the action of external forces perturbing the system during a finite period of time.

A fixed point or stationary point is a state

$$x^{*}=\{x_n^{*}:n=1,2,\dots ,N\},$$

(9)

for which there exists a time $t^{*}$, such that

$$x(t,x^{*})=x^{*}(t\ge t^{*}).$$

(10)

The norm of a dynamic state is denoted as

$$\Vert x(t,x_0)\Vert =\sqrt{\sum _{n=1}^N{x}_n^2(t,x_0)},$$

(11)

which is termed a vector norm or Euclidean norm. This norm describes the character of motion that can be singular, or collapsing, and it can be nonsingular, or non-collapsing.

The motion, starting at a point $x_0$, is strongly singular, if there exists such a finite moment of time $t_c=t_c\left(x_0\right)$, called critical, when the trajectory diverges to infinity,

$$\Vert x(t,x_0)\Vert \to \infty (t\to t_c).$$

(12)

The motion is weakly singular, when the trajectory diverges as time goes to infinity:

$$\Vert x(t,x_0)\Vert \to \infty (t\to \infty ).$$

(13)

The motion is nonsingular, if the trajectory is bounded and never diverges.

2.2. Stability of Solutions

Solutions to evolution equations can be stable or unstable [109]. Finding the regimes of stability for the solutions to the equations of motion is one of the main problems in studying the evolution of social systems. There are several types of stability. The most general type of stability is the Lagrange stability.

Lagrange Stability. The motion, starting at a point $x_0$, is Lagrange stable, if it is nonsingular for all times, that is, the trajectory is always bounded,

$$\Vert x(t,x_0)\Vert <\infty (t\ge t_0).$$

(14)

Global Lagrange Stability. The motion is globally Lagrange stable, if it is Lagrange stable for all initial conditions, so that

$$\Vert x(t,x_0)\Vert <\infty (x_0\in {\mathbb{R}}^N,t\ge t_0).$$

(15)

Poincaré Stability. The solution $x(t,x_0)$ is Poincaré stable if, for any $ϵ>0$, there exist ${\delta }_{ϵ}$ and a time $t_{ϵ}$, called recurrence time, such that, if $|x_0-x_0^{\prime }|<{\delta }_{ϵ}$, then

$$\Vert x(t_{ϵ},x_0^{\prime })-x_0\Vert <{ϵ}^{}.$$

(16)

Since any moment of time can be accepted as a new initial point, there can exist an infinite sequence or recurrence times.

Lyapunov Stability. The solution $x(t,x_0)$ is Lyapunov stable if, for any $ϵ>0$, there exists ${\delta }_{ϵ}$ such that, if $\Vert x_0-x_0^{\prime }\Vert <{\delta }_{ϵ}$, then for all $t>0$,

$$\Vert x(t,x_0)-x(t,x_0^{\prime })\Vert <{ϵ}^{}(t>0).$$

(17)

This means that, if two trajectories, at the initial time, start sufficiently close to each other, they remain close to each other for all times $t>0$.

Global Lyapunov Stability. The solution $x(t,x_0)$ is globally Lyapunov stable, if it is Lyapunov stable for any initial condition $x_0\in {\mathbb{R}}^N$.

Asymptotic Lyapunov Stability. The solution $x(t,x_0)$ is asymptotically stable, if it is Lyapunov stable and there exists ${\delta }_{ϵ}>0$ such that if $\parallel x_0-x_0^{\prime }\parallel <{\delta }_{ϵ}$, then

$$\underset{t\to \infty }{lim}\Vert x(t,x_0)-x(t,x_0^{\prime })\Vert =0.$$

(18)

That is, if two asymptotically stable trajectories, at the initial time, start sufficiently close to each other, they coincide for the time tending to infinity.

Global Asymptotic Stability. The solution $x(t,x_0)$ is globally asymptotically stable, if it is asymptotically stable for any initial condition $x_0\in {\mathbb{R}}^N$.

Definitions of stability can be straightforwardly reformulated for fixed points. For the autonomous equations, fixed points $x^{*}$ are given by the condition,

$$\frac{d{x}^{*}}{dt}=f\left(x^{*}\right)=0.$$

(19)

Lyapunov Stable Fixed Point. The fixed point $x^{*}$ is Lyapunov stable, if for any $ϵ>0$ there exists ${\delta }_{ϵ}>0$ such that, if $\parallel x_0-x^{*}\parallel <{\delta }_{ϵ}^{}$, then

$$\Vert x(t,x_0)-x^{*}\Vert <{ϵ}^{}(t\ge 0)$$

(20)

for all t > 0.

Globally Lyapunov Stable Fixed Point. The fixed point $x^{*}$ is globally Lyapunov stable if it is Lyapunov stable for all initial conditions $x_0\in {\mathbb{R}}^N$.

Asymptotically Stable Fixed Point. The fixed point $x^{*}$ is asymptotically stable, if it is Lyapunov stable and there exists $\delta >0$ such that, if $\parallel x_0-x^{*}\parallel <\delta$, then

$$\underset{t\to \infty }{lim}\Vert x(t,x_0)-x^{*}\Vert =0.$$

(21)

Globally Asymptotically Stable Fixed Point. The fixed point $x^{*}$ is globally asymptotically stable if it is asymptotically stable for all initial conditions $x_0\in {\mathbb{R}}^N$.

An asymptotically stable fixed point is called an attractor, since it attracts all trajectories starting in the vicinity of an initial point $x_0$. The maximal region in the phase space, around an initial point $x_0$, from where all trajectories tend to an attractor $x^{*}$, is termed the basin of attraction. The motion is called regular, when it is Lyapunov stable, and the motion is chaotic, when it is Lyapunov unstable.

2.3. Method of Linearization

Lyapunov developed two methods for analysing the stability of motion. The most often used is the linearisation method that is briefly described just below.

The analysis of stability starts with considering a small deviation $\delta x\left(t\right)$ from the evolution of the variable $x\left(t\right)$:

$$x\left(t\right)⟼x\left(t\right)+\delta x\left(t\right),$$

(22)

where the deviation is a vector,

$$\delta x\left(t\right)=\{\delta x_n\left(t\right):n=1,2,\dots ,N\}.$$

(23)

Substituting this into the evolution equation, to keep in mind is that $x\left(t\right)$ is a solution to this equation. Taking into account only the linear deviation, one gets:

$$\frac{d}{dr}\delta x\left(t\right)=\widehat{J}(x,t)\delta x\left(t\right),$$

(24)

with the Jacobian matrix,

$$\widehat{J}(x,t)=\left[J_{mn}(x,t)\right],J_{mn}\equiv \frac{\partial f_m(x,t)}{\partial x_n}.$$

(25)

One can define the eigenvalues and eigenvectors of the Jacobian matrix by the eigenproblem,

$$\widehat{J}(x,t){\phi }_n(x,t)=J_n(x,t){\phi }_n(x,t),$$

(26)

in which the eigenvectors are the vector-columns,

$${\phi }_n(x,t)=\{{\phi }_{nm}(x,t):m=1,2,\dots ,N\}.$$

(27)

Explicitly, Equation (26) can be written in the form,

$$\sum _m(J_{nm}-J_n{\delta }_{nm}){\phi }_{nm}=0,$$

where, for brevity, the variables are not shown and ${\delta }_{nm}$ is the Kronecker delta. This equation possesses a nontrivial solution, provided that the eigenvalues are given by the condition of the zero determinant,

$$|\widehat{J}(x,t)-J_n(x,t)|=0.$$

(28)

The linearisation is generally performed in the vicinity of fixed points, provided these exist, for which it is required that the evolution equation in the limit of large time becomes autonomous,

$$f(x,t)\simeq f\left(x\right)(t\to \infty ).$$

(29)

Then, the Jacobian, as well as its eigenvalues and eigenvectors, in the limit of large time, do not depend on time:

$$\widehat{J}(x,t)\simeq \widehat{J}\left(x\right),{\phi }_n(x,t)\simeq {\phi }_n\left(x\right)(t\to \infty ).$$

(30)

Considering small deviations,

$$\delta x\left(t\right)=x\left(t\right)-x^{*},$$

(31)

from a fixed point, given by the definition,

$$f\left(x^{*}\right)=0,$$

(32)

one has: the equation,

$$\frac{d}{dt}\delta x\left(t\right)=\widehat{J}\left(x^{*}\right)\delta x\left(t\right).$$

(33)

The Jacobian eigenproblem takes the form,

$$\widehat{J}\left(x^{*}\right){\phi }_n\left(x^{*}\right)=J_n\left(x^{*}\right){\phi }_n\left(x^{*}\right).$$

(34)

The solution to the linear equation for the deviation can be represented as an expansion over the eigenvectors of the Jacobian matrix:

$$\delta x\left(t\right)=\sum _n{c}_{n}exp\left\{J_n\left(x^{*}\right)t\right\}{\phi }_n\left(x^{*}\right).$$

(35)

The Lyapunov exponents, or characteristic exponents, related to a fixed point, are the real parts of the Jacobian eigenvalues evaluated at this fixed point,

$${\lambda }_n\equiv \mathrm{Re}{J}_n\left(x^{*}\right).$$

(36)

The collection $\left\{{\lambda }_n\right\}$ of all Lyapunov exponents forms the Lyapunov spectrum. The number of Lyapunov exponents equals the dimension of the phase space. The largest Lyapunov exponent from the Lyapunov spectrum is called the convergence rate,

$$\lambda \equiv \underset{n}{sup}{\lambda }_n.$$

(37)

Lyapunov exponents characterise the asymptotic stability of fixed points and describe the type of the asymptotic dynamics. The sign of the largest Lyapunov exponent (37) determines the asymptotic convergence or divergence of trajectories with infinitesimally close initial conditions. If the largest exponent is negative, the motion is regular, so that the trajectories with infinitesimally close initial conditions converge to each other. However, if the largest exponent is positive, then the initially close trajectories diverge, and the dynamical system is chaotic. More information on chaotic motion can be found in [110,111,112].

Lyapunov theorem. If all Lyapunov exponents for a fixed point are negative, this fixed point is asymptotically stable. And if at least one of the Lyapunov exponents is positive, this fixed point is asymptotically unstable.

As a simplest case, let us consider a one-dimensional autonomous dynamical system. The Jacobian is just the derivative,

$$J\left(x\right)=\frac{\partial f\left(x\right)}{\partial x}.$$

For each fixed point, there is a single Lyapunov exponent

$$\lambda =\mathrm{Re}J\left(x^{*}\right).$$

The fixed point is asymptotically stable, if $\lambda <0$. Then, there exists a basin of attraction from which the trajectories tend to the given fixed point. The fixed point is termed neutrally stable, when $\lambda =0$. As so, the fixed point can be a centre, around which the trajectory oscillates. The fixed point is asymptotically unstable, when $\lambda >0$. Then, there is no basin of attraction for this fixed point, but all trajectories diverge from it.

2.4. Plane Motion

The motion on a plane is described by a two-dimensional dynamical system. That is, the phase space is ${\mathbb{R}}^2$. For an autonomous system, there are two differential equations:

$$\frac{dx}{dt}=f_1(x,y),\frac{dy}{dt}=f_2(x,y),$$

(38)

with initial conditions,

$$x\left(0\right)=x_0,y\left(0\right)=y_0.$$

(39)

Dividing the second equation over the first gives the relation,

$$\frac{dy}{dx}=\frac{f_2(x,y)}{f_1(x,y)},$$

(40)

defining the trajectory,

$$y=y(x,x_0,y_0).$$

(41)

The set of the trajectories in the plane for different initial conditions forms the phase portrait of the dynamical system,

$$\{y(x,x_0,y_0):x_0\in {\mathbb{R}}^2,y_0\in {\mathbb{R}}^2\}.$$

(42)

The Jacobian is a two-by-two matrix with the elements,

$$\begin{gathered} J_{11}=\frac{\partial f_1}{\partial x},J_{12}=\frac{\partial f_1}{\partial y}, \\ {J}_{21}=\frac{\partial f_2}{\partial x},J_{22}=\frac{\partial f_2}{\partial y}, \end{gathered}$$

(43)

where $J_{mn}\equiv J_{mn}(x,y)$, but, for brevity, the dependence on x and y is omitted. The Jacobian eigenvalues are given by the equation,

$$\begin{array}{c}\hfill \left|\begin{array}{cc}{J}_{11}-J\hfill & J_{12}\hfill \\ J_{21}\hfill & J_{22}-J\hfill \end{array}\right|=0,\end{array}$$

(44)

resulting in two solutions:

$$J_{1,2}=\frac{1}{2}\left[\mathrm{Tr}\widehat{J}\pm \sqrt{{\left(\mathrm{Tr}\widehat{J}\right)}^2-4\det \widehat{J}}\right],$$

(45)

in which

$$\mathrm{Tr}\widehat{J}=J_{11}+J_{22},\det \widehat{J}=J_{11}{J}_{22}-J_{12}{J}_{21}.$$

The solutions for the fixed points, given by the equations,

$$f_1(x^{*},y^{*})=0,f_2(x^{*},y^{*})=0,$$

(46)

can be of the following types depending on the Jacobian eigenvalues $J_n=J_n(x^{*},y^{*})$.

Stable node: $J_1\le J_2<0$, hence ${\lambda }_1\le {\lambda }_2<0$. There exists a basin of attraction, from which all trajectories tend to this fixed point without oscillations.

Unstable node: $J_1\ge J_2>0$, hence ${\lambda }_1\ge {\lambda }_2>0$. There is no basin of attraction, and all trajectories, starting in the vicinity of this fixed point diverge without oscillations.

Stable focus: $J_1=\lambda +i\omega$ and $J_2=\lambda -i\omega$, with $\lambda <0$ and $\omega >0$. There is a basin of attraction, from which all trajectories tend to the fixed point, oscillating around it.

Unstable focus: $J_1=\lambda +i\omega$ and $J_2=\lambda -i\omega$, with $\lambda >0$ and $\omega >0$. There is no basin of attraction, and all trajectories, starting in the vicinity of the fixed point, diverge from the point, oscillating around it.

Elliptic point: $J_1=i\omega$ and $J_2=-i\omega$, with $\omega >0$. There exists a limit cycle, either stable or unstable. A stable limit cycle is a closed orbit, surrounding the fixed point, and attracting all trajectories starting from the related basin of attraction. An unstable limit cycle is a closed orbit, surrounding the fixed point, and repelling all trajectories, starting in its vicinity.

Saddle point: $J_1<0$ and $J_2>0$, hence ${\lambda }_1<0<{\lambda }_2$. There are only two trajectories tending to the point, while all other trajectories diverge from the point.

For continuous dynamical systems on the plane, there exists the Poincare–Bendixson theorem showing that such systems enjoy rather regular motion and do not exhibit chaotic behaviour.

Poincare–Bendixson theorem. If a trajectory of a continuous two-dimensional dynamical system is bounded, that is, Lagrange stable, then it approaches either a fixed point or a limit cycle.

2.5. Evolution Equations

Evolution equations are typically represented by differential equations and sometimes by delay differential equations. The general structure of such equations is mentioned in Sections below and then these equations are specified by several examples.

2.5.1. Differential Rate Equations

Social systems are composed of population groups whose members interact with each other cooperating or competing. The number of members in each group can vary, which can be characterised by differential equations.

Let us consider a society composed of several different groups of populations enumerated by the index $i=1,2,\dots$. In an ith group there are $N_i$ members. The number of the members, $N_i=N_i\left(t\right)$, is not fixed but can vary, e.g., because of the members births and deaths, as well as due to newcomers joining the society from outside. The total number of the society members is

$$N=\sum _i{N}_i=N\left(t\right).$$

(47)

In other cases, N can represent a company capitalisation, with $N_i$ being the capitalisation of the company parts.

The number of the society members can be quite large making it more convenient to use reduced numbers, for example normalising $N_i$ by a fixed scaling $N_0$:

$$x_i\equiv \frac{N_i}{N_0}(i=1,2,\dots ).$$

(48)

Often, $N_0$ can be chosen as the total number of the society members at the initial time, $N\left(0\right)$, however this is not compulsory. The scaling number is chosen according to the convenience and can be different for different cases. The set, $x_i:i=1,2,\dots$, of all variables forms the society state.

The evolution equations for the society groups generally have the form of the rate equation,

$$\frac{d{N}_i}{dt}=R_i{N}_i+{\mathsf{\Phi }}_i,$$

(49)

where $R_i=R_i(x,t)$ is the effective rate of the population group change and $\mathsf{\Phi }=\mathsf{\Phi }(x,t)$ is an external population flux. When the rate $R_i$ is positive, the group population increases, while if the rate $R_i$ is negative, the group population decreases.

The effective rate is often taken in the form of

$$R_i={\gamma }_i+\sum _j{A}_{ij}{N}_j,$$

(50)

where ${\gamma }_i$, e.g., is a birth (death) rate and the second term reflects the influence on the rate of other groups. The equation for the number of the group members reads:

$$\frac{d{N}_i}{dt}=\left({\gamma }_i+\sum _j{A}_{ij}{N}_j\right)N_i+{\mathsf{\Phi }}_i,$$

(51)

and in the reduced form the evolution equation becomes

$$\frac{d{x}_i}{dt}=\left({\gamma }_i+\sum _j{a}_{ij}{x}_j\right)x_i+{\phi }_i,$$

(52)

with the notation,

$$a_{ij}\equiv A_{ij}{N}_0,{\phi }_i\equiv \frac{{\mathsf{\Phi }}_i}{N_0}.$$

2.5.2. Delay Differential Equations

Sometimes, modelling the effective rates of evolution equations, one takes into account that the description of realistic processes involves delayed actions, so that the evolution equation for a species $x_i$ becomes delay-differential equation [113,114]. The delay differential equations belong to the class of functional equations, similarly to partial differential equations, which are infinite dimensional [115].

In the reduced form, the delay equation for the vector,

$$x\left(t\right)=\{x_n\left(t\right):n=1,2,\dots \},$$

with different discrete delays, ${\tau }_n$, for each population, reads:

$$\frac{d}{dt}x\left(t\right)=f(x\left(t\right),x_1(t-{\tau }_1),x_2(t-{\tau }_2),\dots ,t),$$

(53)

where $f=\left\{f_n\right\}$ is also a vector. The history conditions for negative times are to be given by prescribed history functions,

$$x_n\left(t\right)=h_n\left(t\right)(t\le 0),$$

(54)

with $h_n\left(t\right)$ being explicit functions of time.

If for each species, ${\tau }_n=\tau$ is the same, the evolution equation takes the form

$$\frac{d}{dt}x\left(t\right)=f(x\left(t\right),x(t-\tau ),t).$$

(55)

The solution of a delay equation can be calculated following the step-by-step method. For instance, let us consider the case of a single population described by Equation (55), with the history condition, $x\left(t\right)=h\left(t\right)$ for $t\le 0$. One considers the solution separately in the subsequent time intervals $[0,\tau ]$, then $[\tau ,2\tau ]$, then $[2\tau ,3\tau ]$, etc. At the first step, the solution,

$$x\left(t\right)=x^{\left(1\right)}\left(t\right)(0\le t\le \tau ),$$

(56)

is defined by the equation,

$$\frac{d}{dt}{x}^{\left(1\right)}\left(t\right)=f(x^{\left(1\right)}\left(t\right),h\left(t\right),t),$$

(57)

with the initial condition,

$$x^{\left(1\right)}\left(0\right)=h\left(0\right).$$

(58)

At the second step, denoting

$$x\left(t\right)=x^{\left(2\right)}\left(t\right)(\tau \le t\le 2\tau ),$$

(59)

one finds the solution from the equation,

$$\frac{d}{dt}{x}^{\left(2\right)}\left(t\right)=f(x^{\left(2\right)}\left(t\right),x^{\left(1\right)}\left(t\right),t)$$

(60)

with the initial condition,

$$x^{\left(2\right)}\left(\tau \right)=x^{\left(1\right)}\left(\tau \right).$$

(61)

Continuing this procedure for subsequent time intervals, at the step k, the solution,

$$x\left(t\right)=x^{\left(k\right)}\left(t\right)((k-1)\tau \le t\le k\tau ),$$

(62)

is given by the equation,

$$\frac{d}{dt}{x}^{\left(k\right)}\left(t\right)=f(x^{\left(k\right)}\left(t\right),x^{(k-1)}\left(t\right),t),$$

(63)

with the initial condition,

$$x^{\left(k\right)}\left((k-1)\tau \right)=x^{(k-1)}\left((k-1)\tau \right).$$

(64)

In practice, the solution at each step is calculated numerically.

2.6. Examples of Evolution Equations

There exist many variants of evolution equations. Several most known equations are briefly dicussed just below.

2.6.1. Malthus Equation

Malthus [116] was interested in the Earth’s population growth. He noticed that the Earth’s population varies according to the law,

$$\frac{dN}{dt}=\gamma N.$$

(65)

The variation rate, $\gamma$, is the difference between the growth, ${\gamma }_{+}$, and death, ${\gamma }_{-}$, rates of the population,

$$\gamma ={\gamma }_{+}-{\gamma }_{-}.$$

(66)

Taking for the scaling $N_0=N\left(0\right)$, one gets: the equation in the reduced form,

$$\frac{dx}{dt}=\gamma x,$$

(67)

with the initial condition,

$$x\left(0\right)=1.$$

(68)

This gives the evolution law,

$$x=e^{\gamma t}.$$

(69)

As is understandable, with time, the population either diminishes to zero, if $\gamma <0$, or tends to infinity, if $\gamma >0$. Since, for the people on the Earth, the birth rate surpasses the death rate, hence $\gamma >0$, the population will grow to infinity. The singular behaviour of the population results in overpopulation, which is called the Malthusian catastrophe. The population explosion leads to wars, epidemics, and hunger. To avoid overpopulation, Malthus advised to keep $\gamma$ equal to zero.

2.6.2. Neo-Malthusian Catastrophe

Some authors have suggested that the reality is even more complicated than that predicted by Malthus, and the Earth population can become strongly singular, diverging at a finite moment of time. This scenario can be compared with the proliferation of cancer cells killing an organism at finite time [117].

The finite-time divergence can happen when, in addition to the birth rate, one includes into the effective rate the influence of population cooperation, so that the rate becomes

$$R=\gamma +AN(\gamma >0,A>0).$$

(70)

Then, the evolution equation reads:

$$\frac{dN}{dt}=(\gamma +AN)N.$$

(71)

With the scale $N_0=N\left(0\right)$, one has the reduced equation,

$$\frac{dx}{dt}=(\gamma +ax)x.$$

(72)

The solution,

$$x=\frac{\gamma e^{\gamma t}}{\gamma -a(e^{\gamma t}-1)}$$

(73)

shows that the population diverges as

$$x\simeq \frac{1}{a(t_c-t)}(t\to t_c-0),$$

(74)

at the finite time,

$$t_c=\frac{1}{\gamma }ln\left(1+\frac{\gamma }{a}\right).$$

(75)

For the Earth’s population, this time corresponds to the year 2027. Thus, the human population is analogous to cancer cells destroying their habitat.

2.6.3. Logistic Equation

Verhulst [118,119] suggested that no population can grow indefinitely, but there exists a carrying capacity, K, above which the population growth is impossible. He proposed the logistic equation,

$$\frac{dN}{dt}=\gamma \left(1-\frac{N}{K}\right)N(\gamma >0,K>0),$$

(76)

taking into account the existence of the carrying capacity K entering the term characterising the competition between the society members. The Malthus equation can be valid only at the very beginning of the population growth when $N\ll K$.

The logistic equation has found a number of applications describing the growth of population in demography, company capitalisation in economics, neuron signals in neural networks, tumuor size in medicine, reactant mass in chemistry, etc. The reduced form of the equation, with the population normalised to $N_0=N\left(0\right)$, is

$$\frac{dx}{dt}=(\gamma -ax)x,$$

(77)

where

$$a\equiv \gamma \frac{N_0}{K}>0.$$

The solution to the logistic equation is the logistic function,

$$x=\frac{\gamma e^{\gamma t}}{\gamma +a(e^{\gamma t}-1)},$$

(78)

also called the sigmoid function. Its limit with respect to time is finite,

$$\underset{t\to \infty }{lim}x=\frac{\gamma }{a}=\frac{K}{N_0}(\gamma >0),$$

(79)

hence there is no singular behaviour.

From the point of view of the Lyapunov analysis, fixed points are defined by the equation,

$$(\gamma -a{x}^{*})x^{*}=0,$$

resulting in two fixed points, $x_1^{*}=0$ and $x_2^{*}=\gamma /a$. The Jacobian,

$$J\left(x\right)=\gamma -2ax,$$

takes the values,

$$J\left(x_1^{*}\right)=\gamma ,J\left(x_2^{*}\right)=-\gamma .$$

Hence, the trivial fixed point, $x_1^{*}=0$, is not stable, but the nontrivial fixed point, $x_2^{*}>0$, is stable.

2.6.4. Lotka-Volterra Model

This model [120,121] describes the coexistence of two populations, where one of them corresponds to prey (x) and the other, to predators (y), because of which it is also named the predator-prey model. This model is also employed in many other applications. For example, in combustion theory these can be a passive radical (x) and an active radical (y), in medicine, it can be susceptible to infection individuals (x) and infective individuals (y), in economics, buying population (x) and sold goods (y).

In the reduced form, the corresponding equations are

$$\frac{dx}{dt}={\gamma }_{1}x-a_{12}yx,\frac{dy}{dt}=-{\gamma }_{2}y+a_{21}xy,$$

(80)

where all parameters are assumed to be positive. The initial conditions are denoted as

$$x\left(0\right)=x_0,y\left(0\right)=y_0.$$

(81)

When $x_0=0$, then $y\left(t\right)$ dies down to zero as time increases. If $y_0=0$, then $x\left(t\right)$ explodes to infinity with time. More interesting are the nonzero initial conditions.

Let us employ the Lyapunov stability analysis, as described in Section 2.3. The fixed-point equations,

$$({\gamma }_1-a_{12}{y}^{*})x^{*}=0,(-{\gamma }_2+a_{21}{x}^{*})y^{*}=0,$$

(82)

yield two fixed points: one is trivial,

$$x_1^{*}=0,y_1^{*}=0,$$

(83)

and the other is nontrivial,

$$x_2^{*}=\frac{{\gamma }_2}{a_{21}},y_2^{*}=\frac{{\gamma }_1}{a_{12}}.$$

(84)

The Jacobian matrix is composed of the elements,

$$J_{11}={\gamma }_1-a_{12}y,J_{12}=-a_{12}x,J_{21}=a_{21}y,J_{22}=-{\gamma }_2+a_{21}x,$$

(85)

which lead to the Jacobian trace and determinant,

$$\mathrm{Tr}\widehat{J}(x,y)={\gamma }_1-{\gamma }_2+a_{21}x-a_{12}y,$$

$$\det \widehat{J}(x,y)=({\gamma }_1-a_{12}y)(-{\gamma }_2+a_{21}y)+a_{12}{a}_{21}xy.$$

At the trivial fixed point, one has:

$$\mathrm{Tr}\widehat{J}(x_1^{*},y_1^{*})={\gamma }_1-{\gamma }_2,\det \widehat{J}(x_1^{*},y_1^{*})=-{\gamma }_1{\gamma }_2.$$

The Jacobian eigenvalues are

$$J_1(x_1^{*},y_1^{*})={\gamma }_1,J_2(x_1^{*},y_1^{*})=-{\gamma }_2.$$

(86)

Hence, the trivial fixed point is unstable, being a saddle point.

For the second fixed point, one gets:

$$\mathrm{Tr}\widehat{J}(x_2^{*},y_2^{*})=0,\det \widehat{J}(x_2^{*},y_2^{*})={\gamma }_1{\gamma }_2,$$

hence, the Jacobian eigenvalues become

$$J_1(x_2^{*},y_2^{*})=i\omega ,J_2(x_2^{*},y_2^{*})=-i\omega ,$$

(87)

where

$$\omega =\sqrt{{\gamma }_1{\gamma }_2}.$$

(88)

This means that the nontrivial fixed point is an elliptic point for any nonzero parameters of the evolution equation. The corresponding temporal behaviour is given by permanent oscillations. The predators and prey oscillate with the same frequency $\omega$, the predator curve being shifted to the right with respect to the prey curve.

There are several other types of equations describing population evolution, some of which are mentioned below.

2.6.5. Singular Malthus Equation

The Malthus Equation (65) can be modified [122] to the form,

$$\frac{dN}{dt}=\gamma N^m(m\ge 1),$$

(89)

whose solution diverges by power law,

$$N=\frac{1}{{\left({ϵ}^{\gamma }\right)}^{1/ϵ}{(t_c-t)}^{1/ϵ}},$$

(90)

if $\epsilon =m-1>0$ and the critical time is

$$t_c=\frac{1}{N{\left(0\right)}^{ϵ}ϵ\gamma }.$$

(91)

For $m=1$, one returns to the original Malthus Equation (65) that yields an exponentially divergent solution at $t\to \infty$.

Such strongly singular solutions were applied to rationalise the super-exponential growth of the human world population, financial markets, material failures, earthquakes, climate changes, and dynamics of other systems (see [123]).

2.6.6. Generalized Lotka-Volterra Model

The Lotka–Volterra model (80) can also be generalised to the case of multiple species for the populations $N_i$, with $i=1,2,\dots$, yielding:

$$\frac{d{N}_i}{dt}=\left({\gamma }_i+\sum _j{A}_{ij}{N}_j\right)N_i.$$

(92)

The dimensionless populations can be introduced by normalising $N_i$, for example, to the total initial population $N\left(0\right)={\sum }_i{N}_i\left(0\right)$, which gives the reduced form,

$$\frac{d{x}_i}{dt}=\left({\gamma }_i+\sum _j{a}_{ij}{x}_j\right)x_i,$$

(93)

in which

$$x_i\equiv \frac{N_i}{N\left(0\right)},a_{ij}\equiv A_{ij}N\left(0\right).$$

The signs of the parameters ${\gamma }_i$ and $a_{ij}$ can be different, providing a large variety of possible solutions [124]. For two species, one gets the original Lotka–Volterra model (80).

2.6.7. Predator–Prey Kolmogorov Model

A particular form generalising the Lotka–Volterra Model, in which

$$\frac{d{N}_1}{dt}=f_1(N_1,N_2)N_1,\frac{d{N}_2}{dt}=f_2(N_1,N_2)N_2,$$

(94)

under the conditions

$$\frac{\partial f_1}{\partial N_2}>0,\frac{\partial f_2}{\partial N_1}<0,$$

(95)

is called the predator–prey Kolmogorov model [125,126]. This model can again be rewritten in the reduced form

$$\frac{d{x}_1}{dt}=f_1(x_1,x_2)x_1,\frac{d{x}_2}{dt}=f_2(x_1,x_2)x_2.$$

(96)

To model concrete cases, it is required to specify the functions $f_1$ and $f_2$.

2.6.8. Jacob–Monod Equations

The Jacob–Monod equations describe a single type of population $N_1$ that is being fed on the nutrient of amount $N_2$,

$$\frac{d{N}_1}{dt}=f\left(N_2\right)N_1,\frac{d{N}_2}{dt}=-\gamma f\left(N_2\right)N_1.$$

(97)

For instance, this can be bacteria, playing the role of predators that are fed on a nutrient playing the role of the prey. The nutrient is getting depleted being consumed by the predators. At the same time, the nutrient is supplied into the system from outside, which is described by the function $f\left(N_2\right)$. The supply function is taken in the form, such that the nutrient becomes depleted ($N_2\to 0$), as time increases, while the bacteria population reaches a fixed point characterising a stationary value [127].

2.6.9. Hutchinson Delayed Equation

There are as well several generalisations of the logistic Equation (76). Thus Ref. [128] considered an effective reproductive rate delayed in time, which gives the equation,

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

(98)

Here K is a fixed carrying capacity and $\tau$ is a delay time. The solution to this equation displays oscillations above the logistic curve.

2.6.10. Multiple-Processes Delayed Equation

There exist many variants of the delay logistic equations [113,129] designed for the description of single-species population N. For example, one can introduce multiple carrying capacities, $K_j$, and multiple delay times, ${\tau }_j$, characterising different processes in the dynamics of a population. The single-species population dynamics, with multiple processes reads:

$$\frac{dN\left(t\right)}{dt}=\gamma \left[1-\sum _j\frac{N(t-{\tau }_j)}{K_j}\right]N\left(t\right).$$

(99)

2.6.11. Peschel-Mende Hyperlogistic Equation

The accelerated growth of population, such that exists for the human world population, can be described [130] by means of two additional powers, m and n, which leads to the equation,

$$\frac{dN}{dt}=\gamma {\left(1-\frac{N}{K}\right)}^n{N}^m.$$

(100)

This equation leads to a solution that can be fitted well to the world population dynamics for some finite intervals of time. The solution to this equation is similar to a modified sigmoid curve, which does not exceeds the given carrying capacity K. The Verhulst logistic Equation (76) is recovered for $m=1$ and $n=1$.

2.6.12. Hyperlogistic Delayed Equations

The Peschel-Mende hyperlogistic equation can be extended to the delayed equation [131],

$$\frac{dN\left(t\right)}{dt}=\gamma {\left[1-\frac{N(t-\tau )}{K}\right]}^n{N}^m\left(t\right).$$

(101)

This equation can also be treated as an extension of the Hutchinson delayed Equation (98). The solution to this delayed equation describes a population that can exceed the fixed carrying capacity K for some finite period of time, although finally it decreases below K.

2.7. Replicator Equation

Quite a popular equation employed for modelling the evolution of different types of social and biological systems is the replicator equation describing population dynamics, including dynamics in genetic theory and in evolution game theory. The replicator equation characterises the ensemble of populations, $N_i$, of different species, each being endowed with a fitness $f_i$. It is possible to consider either discrete or continuous in time processes. The continuous replicator equation is

$$\frac{d{N}_i}{dt}=(f_i-\overline{f})N_i,$$

(102)

where

$$\overline{f}\equiv \frac{1}{N}\sum _i{f}_i{N}_i$$

(103)

is the average society fitness. The total population,

$$N=\sum _i{N}_i=\mathrm{const},$$

(104)

is assumed to be constant, hence the equation is defined on a simplex [124].

In terms of the reduced fractions,

$$x_i\equiv \frac{N_i}{N}\left(\sum _i{x}_i=1\right),$$

(105)

the equation takes the form

$$\frac{d{x}_i}{dt}=(f_i-\overline{f})x_i,$$

(106)

with the average fitness

$$\overline{f}\equiv \sum _i{f}_i{x}_i.$$

(107)

This equation leads to a conclusion that in a society of several species only the species with a largest fitness survives.

To illustrate how the equation works, let us consider quite a simple case of two populations, $N=2$, so that $x_1+x_2=1$. Then the system of two equations can be reduced to a single equation, e.g.,

$$\frac{d{x}_1}{dt}=(f_1-f_2)(1-x_1)x_1.$$

(108)

There are two fixed points, $x_1^{*}=0$ and $x_1^{*}=1$. The Jacobian is

$$J\left(x_1\right)=(f_1-f_2)(1-2x_1).$$

(109)

At the fixed points, this gives

$$J\left(0\right)=f_1-f_2,J\left(1\right)=-(f_1-f_2).$$

(110)

Therefore, if the fitness of the first species $f_1$ is larger than that of the second species, $f_2$, then the stable fixed point corresponds to

$$x_1^{*}=1(f_1>f_2).$$

(111)

However, when the fitnesses are such that $f_1$ is smaller than $f_2$, then only the second species survives, since

$$x_1^{*}=0(f_1<f_2).$$

(112)

This result one often formulates as a sentence “the strongest survives”.

2.8. Free Replicator Equation

It is interesting that the same qualitative result follows from the replicator equation that is free from the normalisation conditions (104) or (105), which can be called free replicator equation. Then one needs to deal with the number of equations equal to the number of species.

Thus, for the case of two species, one has two equations,

$$\frac{d{x}_1}{dt}=(f_1-x_1{f}_1-x_2{f}_2)x_1,\frac{d{x}_2}{dt}=(f_2-x_1{f}_1-x_2{f}_2)x_2.$$

(113)

Let us assume that the species fitnesses are different, $f_1\ne f_2$, otherwise the species could not be distinguished.

There are three fixed points, where either

$$x_1^{*}=0,x_2^{*}=0,$$

(114)

or

$$x_1^{*}=0,x_2^{*}=1,$$

(115)

or

$$x_1^{*}=1,x_2^{*}=0.$$

(116)

The Jacobian matrix has the components,

$$\begin{gathered} J_{11}(x_1,x_2)=(1-2x_1)f_1-x_2{f}_2,J_{12}(x_1,x_2)=-x_1{f}_2, \\ {J}_{21}(x_1,x_2)=-x_2{f}_1,J_{22}(x_1,x_2)=(1-2x_2)f_2-x_1{f}_1. \end{gathered}$$

(117)

At the fixed point (114), one has:

$$J_{11}(0,0)=f_1,J_{12}(0,0)=0,J_{21}(0,0)=0,J_{22}(0,0)=f_2.$$

(118)

This gives the positive Lyapunov exponents,

$${\lambda }_1=f_1,{\lambda }_2=f_2,$$

(119)

showing that this fixed point is not stable.

At the fixed point (115), the Jacobian components are

$$J_{11}(0,1)=f_1-f_2,J_{12}(0,1)=0,J_{21}(0,1)=-f_1,J_{22}(0,1)=-f_2,$$

(120)

which gives the Lyapunov exponents,

$${\lambda }_1=f_1-f_2,{\lambda }_2=-f_2.$$

(121)

Hence the fixed point (115) is stable when $f_1<f_2$. Then only the second species survives.

At the fixed point (116), for the Jacobian one gets:

$$J_{11}(1,0)=-f_1,J_{12}(1,0)=-f_2,J_{21}(1,0)=0,J_{22}(1,0)=f_2-f_1,$$

(122)

thence the Lyapunov exponents are

$${\lambda }_1=f_2-f_1,{\lambda }_2=-f_1.$$

(123)

This tells us that the fixed point (116) is stable provided that $f_1>f_2$, when the second species dies out, while the first species survives.

2.9. Influence of Noise

The evolution equations considered above are called deterministic, since, given the form of an equation and initial conditions, the following dynamics is uniquely defined. In reality, it may happen that the evolution of a society is subject to random perturbations, termed random noise. Then, instead of Equation (7), one comes to the equation,

$$\frac{dx}{dt}=f(x,t)+g(x,t)\xi \left(t\right),$$

(124)

containing, in addition to $f(x,t)$, a term describing the action of an external noise. Generally, setting

$$\xi \left(t\right)dt=dW\left(t\right),$$

(125)

one writes Equation (124) in the form

$$dx=f(x,t)dx+g(x,t)dW\left(t\right),$$

(126)

since the noise term is often modelled by the Gaussian white noise that is not differentiable. In the case of white noise, the variable $W\left(t\right)$ characterises a standard Wiener process. For time in the interval $[0,t]$, the random variable $W\left(t\right)$ is drawn from the normal law with zero mean and the standard deviation $\sqrt{t}$. Equation (126) is called stochastic differential equation (see Refs. [132,133,134,135]). The quantity $f(x,t)$ is called drift term and $g(x,t)$, diffusion term. For a while, a single variable $x=x\left(t\right)$ is considered. The generalisation to many variables is given at the end of this Section.

In practical calculations, the dynamics of $x\left(t\right)$, in the presence of noise, is often represented through the Euler–Maruyama scheme [136]. For this purpose, the interval $[0,t]$ is partitioned in discrete points $t_n=n\Delta t$, with $n=0,1,2,\dots$. Then Equation (126) is integrated between $t_n$ and $t_{n+1}=t_n+\Delta t$, which gives

$$x\left(t_{n+1}\right)-x\left(t_n\right)={\int }_{t_n}^{t_n+\Delta t}f(x,t)dt+{\int }_{t_n}^{t_n+\Delta t}g(x,t)dW\left(t\right).$$

(127)

The first integral in the right-hand side is a known Riemann integral that, for small $\Delta t$, can be written as

$${\int }_{t_n}^{t_n+\Delta t}f(x,t)dt\cong f(x\left(t_n\right),t_n)\Delta t.$$

(128)

The Riemann integral does not depend on the choice of the variable t inside the interval $[t_n,t_n+\Delta t]$, so that instead of $t_n$ in the right-hand side of the above formula it is possible to take any $t_n^{\prime }\in [t_n,t_n+\Delta t]$. However, the integral over a stochastic variable depends on this choice. There exist two accepted ways of choosing the points of t in the second integral of Equation (127), following Stratonovich or Ito [132,133,134,135]). According to Stratonovich, for infinitesimally small $\Delta t$, one takes $[g(x(t_n+\Delta ),t_n+\Delta )+g(x\left(t_n\right),t_n)]/2$, while according to Ito one has to take $g(x\left(t_n\right),t_n)$. When employing the Euler-Maruyama numerical scheme [136], one uses the Ito representation. Then the second integral in Equation (127), understood in the sense of Ito, then called the Ito integral, is written as

$${\int }_{t_n}^{t_n+\Delta t}g(x,t)dW\left(t\right)=g(x\left(t_n\right),t_n)\left[W(t_n+\Delta t)-W\left(t_n\right)\right].$$

(129)

The random variable, $W\left(t_n\right)$, is drawn from the normal distribution with zero mean and standard deviation $\sqrt{\Delta t}$. For each $t_1<t_2$, the normal random variable, $W\left(t_2\right)-W\left(t_1\right)$, is independent of the random variable $W\left(t_1\right)$. The difference can be represented as

$$W(t_n+\Delta t)-W\left(t_n\right)=R\left(t_n\right)\sqrt{\Delta t},$$

(130)

where $R\left(t_n\right)$ is drawn from the normal distribution with zero mean and standard deviation one. Thus, one comes to the iterative equation,

$$x\left(t_{n+1}\right)=x\left(t_n\right)+f(x\left(t_n\right),t_n)\Delta t+g(x\left(t_n\right),t_n)R\left(t_n\right)\sqrt{\Delta t}.$$

(131)

In those cases, where the diffusion term is constant, $g(x,t)=\sigma$, one has:

$$x\left(t_{n+1}\right)=x\left(t_n\right)+f(x\left(t_n\right),t_n)\Delta t+\sigma R\left(t_n\right)\sqrt{\Delta t}.$$

(132)

With several variables $x_i$, where $i=1,2,\dots$, and several noise terms, one has:

$$d{x}_i=f_i(x_i,t)dt+\sum _j{g}_{ij}(x,t)d{W}_j\left(t\right).$$

(133)

For a diagonal diffusion matrix $g_{ij}(x,t)={\sigma }_i{\delta }_{ij}$, this reduces to

$$x_i\left(t_{n+1}\right)=x_i\left(t_n\right)+f_i(x\left(t_n\right),t_n)\Delta t+{\sigma }_{i}R\left(t_n\right)\sqrt{\Delta t}.$$

(134)

The existence of noise superimposes on a dynamical trajectory random deviations, whose amplitude depends on the value of $\sigma$.

2.10. Fokker-Planck Equation

In the presence of noise, the system trajectory is not uniquely defined, but there exists a set of possible trajectories depending on the particular realisation of the random noise. For the stochastic differential Equation (126), the probability density of x is described by the Fokker–Planck equation [137],

$$\frac{\partial }{\partial t}p(x,t)=-\frac{\partial }{\partial x}\left[f(x,t)p(x,t)\right]+\frac{{\partial }^2}{\partial x^2}\left[D(x,t)p(x,t)\right],$$

(135)

with the diffusion coefficient,

$$D(x,t)\equiv \frac{1}{2}{g}^2(x,t).$$

(136)

Introducing the probability current,

$$j(x,t)\equiv f(x,t)p(x,t)-\frac{\partial }{\partial x}\left[D(x,t)p(x,t)\right]$$

(137)

yields the continuity equation,

$$\frac{\partial }{\partial t}p(x,t)+\frac{\partial }{\partial x}j(x,t)=0.$$

(138)

One often is interested in the stationary solution to the Fokker–Planck equation, provided it exists. A stationary solution can arise as time tends to infinity, so that the functions $f(x,t)$, $g(x,t)$, and $D(x,t)$ tend to time-independent expressions $f\left(x\right)$, $g\left(x\right)$, and $D\left(x\right)$. Respectively, the probability density tends to $p\left(x\right)$, satisfying the equation

$$\frac{d}{dx}j\left(x\right)=0,$$

(139)

in which the stationary probability current is

$$j\left(x\right)=f\left(x\right)p\left(x\right)-\frac{d}{dx}\left[D\left(x\right)p\left(x\right)\right].$$

(140)

The latter equation implies that

$$j\left(x\right)=\mathrm{const}.$$

(141)

Suppose that the quantity x varies in an interval starting at $x_0$, so that $x\ge x_0$. Assuming the reflecting boundary condition,

$$j\left(x_0\right)=0,$$

(142)

one gets:

$$j\left(x\right)=0,$$

(143)

which is equivalent to the equation,

$$f\left(x\right)p\left(x\right)-\frac{d}{dx}\left[D\left(x\right)p\left(x\right)\right]=0.$$

(144)

The solution to the latter equation is

$$p\left(x\right)=\frac{C}{D\left(x\right)}exp\left\{{\int }_{x_0}^x\frac{f\left(x^{\prime }\right)}{D\left(x^{\prime }\right)}d{x}^{\prime }\right\},$$

(145)

with the constant C defined from the normalisation condition,

$$\int p\left(x\right)dx=1,$$

(146)

where the integration is over the whole range of variation of x.

Expression (145) is named potential solution, as far as it can be rewritten in the form,

$$p\left(x\right)=Cexp\{-U(x\left)\right\},$$

(147)

with the potential,

$$U\left(x\right)=lnD\left(x\right)-{\int }_{x_0}^x\frac{f\left(x^{\prime }\right)}{D\left(x^{\prime }\right)}d{x}^{\prime }.$$

(148)

All above formulas can be easily generalised to a multivariate case, where the Fokker-Planck equation reads:

$$\frac{\partial }{\partial t}p(x,t)=-\sum _i\frac{\partial }{\partial x_i}\left[f_i(x,t)p(x,t)\right]+\sum _{ij}\frac{{\partial }^2}{\partial x_i\partial x_j}\left[D_{ij}(x,t)p(x,t)\right],$$

(149)

with the diffusion matrix of the elements,

$$D_{ij}(x,t)=\frac{1}{2}\sum _k{g}_{ik}(x,t)g_{jk}(x,t).$$

(150)

3. Generalised Evolution Equations

Often, in the evolution equations the rates of change and population interactions are treated as constant parameters. In more realistic, and hence more general, scenarios, the rates may be functions of the population variables. Several such generalised equations are considered below.

3.1. Functional Carrying Capacity

In addition to direct interactions of different populations, there exist indirect interactions through the influence of the populations on the mutual carrying capacities. This concerns as well the influence of populations on their own carrying capacity.

For example, the carrying capacity of a human society depends on the activity of humans. The technological evolution of humans has allowed increase of effective carrying capacity for humans. At the same time, humans can destroy their carrying capacity, e.g., by destroying their habitat.

To formulate the equations with a functional carrying capacity, let us start with the standard form of an evolution equation,

$$\frac{d{N}_i}{dt}=\left({\gamma }_i+\sum _j{A}_{ij}{N}_j\right)N_i+{\mathsf{\Phi }}_i,$$

(151)

where ${\gamma }_i$ is an effective birth rate for a population, if ${\gamma }_i$ is positive, or a death rate, when ${\gamma }_i$ is negative. If one considers not population evolution, but the dynamics of production of a firm or such, then a positive ${\gamma }_i$ takes the sense of gain rate, while a negative ${\gamma }_i$ means loss rate. The quantity $A_{ij}$ describes effective interactions between the ith and jth populations. Commonly, $A_{ij}$ is treated as a constant or, sometimes, as a given function fixed by external forces [138,139,140].

Generally, in the population interactions it is admissible to distinguish direct interactions and indirect interactions through the mutual influence on their carrying capacities [123,141]. Therefore, the effective interactions can be represented in the form,

$$A_{ij}=-\frac{B_{ij}}{K_j}(K_j>0),$$

(152)

with $B_{ij}$ being a direct interaction parameter and $K_j$, carrying capacity of the jth population. The influence of populations on the carrying capacity is assumed, which implies that it depends on these populations:

$$K_j=K_j(N_1(t-{\tau }_1),N_2(t-{\tau }_2),\dots ).$$

(153)

The populations entering $K_j$ are delayed, since to induce a change in the carrying capacity requires time. The term ${\mathsf{\Phi }}_i$ is caused by the influx of populations from outside, if any.

It is more convenient to deal with reduced dimensionless quantities, for which we introduce

$$x_i\equiv \frac{N_i}{N_0},{\phi }_i\equiv \frac{{\mathsf{\Phi }}_i}{N_0}.$$

(154)

Here, the normalisation quantity $N_0$, for a while, is not defined. It is chosen below so that to simplify the reduced equation. Then the evolution equation reads:

$$\frac{d{x}_i}{dt}=\left({\gamma }_i-\sum _j\frac{B_{ij}{N}_0}{K_j}{x}_j\right)x_i+{\phi }_i.$$

(155)

If there are no external fluxes from outside, the term ${\phi }_i$ should be omitted.

To illustrate the structure of the evolution equation, let us study the case with a single population without external flux. Then Equation (155) reduces to

$$\frac{dx}{dt}=\left(\gamma -\frac{B{N}_0}{K}x\right)x.$$

(156)

Using the sign function,

$$\begin{array}{c}\hfill \mathrm{sgn}x\equiv \left\{\begin{array}{cc}\hfill -1,& x<0,\hfill \\ \hfill 0,& x=0,\hfill \\ \hfill 1,& x>0,\hfill \end{array}\right.\end{array}$$

makes it possible to write:

$$\gamma =\left|\gamma \right|\mathrm{sgn}\gamma ,B=\left|B\right|\mathrm{sgn}B.$$

(157)

Defining the dimensionless carrying capacity,

$$q\equiv \frac{K}{N_0}\left|\frac{\gamma }{B}\right|,$$

(158)

and measuring time in units of ${\gamma }^{-1}$ yields the evolution equation,

$$\frac{dx}{dt}=\left(\mathrm{sgn}\gamma \right)x-\left(\mathrm{sgn}B\right)\frac{x^2}{q}.$$

(159)

Formally, this looks as a logistic equation. However the key difference here is the functional dependence of the carrying capacity on the delayed population,

$$q=q\left(x\right(t-\tau \left)\right).$$

(160)

To find out the explicit form of the functional carrying capacity, let us assume that it can be expanded in powers of the population, so that

$$K\simeq K_0+K_{1}N+K_2{N}^2+\dots .$$

This is equivalent to the expansion,

$$q\left(x\right)\simeq q_0+q_{1}x+q_2{x}^2+\dots .$$

(161)

The carrying capacity has to retain its sense, hence to be finite for asymptotically small amount of population. This means that

$$q_0=\frac{K_0}{N_0}\left|\frac{\gamma }{B}\right|>0.$$

(162)

To simplify the equation, one may choose the normalisation quantity $N_0$ in the form,

$$N_0=K_0\left|\frac{\gamma }{B}\right|;$$

(163)

hence, $q_0=1$, as a result of which the expansion of the carrying capacity becomes

$$q\left(x\right)\simeq 1+b_{1}x+b_2{x}^2+\dots .$$

(164)

When the influence of the population on the carrying capacity is small, it is possible to limit ourselves by the linear approximation, $1+b_{1}x$, as has been considered in Ref. [141]. However, this form, when used under sufficiently strong destructive influence of population on the carrying capacity, can lead to the appearance of zero in the denominator and the occurrence of unreasonable divergencies.

To find a more general expression for the carrying capacity, valid under arbitrarily strong influence of populations, it is necessary to find an effective limit of expansion (161). For this purpose, one can resort to the exponential summation guaranteeing a positive effective limit [142,143,144]. This gives:

$$q\left(x\right)=exp\left\{bx\right(t-\tau \left)\right\}.$$

(165)

In this way, the evolution equation reduces to [145]

$$\frac{dx}{dt}=\left(\mathrm{sgn}\gamma \right)x-\left(\mathrm{sgn}B\right)x^{2}exp\{-bx(t-\tau )\}.$$

(166)

The production parameter, b, characterises the influence of the population on the carrying capacity. The production parameter ($b>0$) is positive for the case of productive activity of population, creating additional means for survival. And the production parameter ($b<0$) becomes negative when the population destroys the given carrying capacity. In such a picture, b is, actually, a destruction parameter.

The initial condition for the delay equation is

$$x\left(t\right)=x_0(t\le 0).$$

Depending on the signs of the parameters $\gamma$ and B, there can happen the following four different types of evolution models characterised by:

$$\mathrm{sgn}\gamma =1,\mathrm{sgn}B=1(\mathrm{gain}+\mathrm{competition}),$$

$$\mathrm{sgn}\gamma =1,\mathrm{sgn}B=-1(\mathrm{gain}+\mathrm{cooperation}),$$

$$\mathrm{sgn}\gamma =-1,\mathrm{sgn}B=1(\mathrm{loss}+\mathrm{competition}),$$

$$\mathrm{sgn}\gamma =-1,\mathrm{sgn}B=-1(\mathrm{loss}+\mathrm{cooperation}).$$

3.2. Evolutionary Stable States

Delay equations allow for the Lyapunov stability analysis, similar to that for the standard differential equations [129]. The fixed points, or stationary states, are defined by the equation,

$$\left(\mathrm{sgn}\gamma \right)x^{*}-\left(\mathrm{sgn}B\right){\left(x^{*}\right)}^{2}exp(-b{x}^{*})=0.$$

(167)

There always exists the trivial solution,

$$x_1^{*}=0(-\infty <b<\infty )$$

(168)

existing for any $\mathrm{sgn}\gamma$ and $\mathrm{sgn}B$.

Nontrivial solutions require the validity of the relation,

$$\frac{\mathrm{sgn}\gamma }{\mathrm{sgn}B}=x^{*}exp(-b{x}^{*})>0,$$

(169)

which imposes the constraint

$$\mathrm{sgn}\gamma =\mathrm{sgn}B(x^{*}>0).$$

(170)

Under this condition, fixed points are given by the equation,

$$x^{*}=exp\left(b{x}^{*}\right).$$

(171)

When the population destroys its carrying capacity, that is $b<0$, there can exist one fixed point in the range,

$$0<x_2^{*}\le 1(b\le 0).$$

(172)

If the population increases its carrying capacity, hence $b>0$, but so that $b<1/e$, there can occur two fixed points, one in the range,

$$1<x_2^{*}<e\left(0<b<\frac{1}{e}\right)$$

(173)

and the other for

$$x_3^{*}>e\left(0<b<\frac{1}{e}\right).$$

(174)

At the bifurcation point $b=1/e$, the fixed points $x_2^{*}$ and $x_3^{*}$ coincide:

$$x_2^{*}=x_3^{*}=e\left(b=\frac{1}{e}\right).$$

(175)

There are no stationary nontrivial states for $b>1/e$.

The fixed-point stability is characterised by the behaviour of small deviations from the fixed point. Substituting into the evolution Equation (166) the definition,

$$x=x^{*}+\delta x,$$

(176)

yields the linearised equation,

$$\frac{d}{dt}\delta x\left(t\right)=\left[\mathrm{sgn}\gamma -2\left(\mathrm{sgn}B\right)x^{*}{e}^{-b{x}^{*}}\right]\delta x\left(t\right)+b\left(\mathrm{sgn}\gamma \right)x^{*}\delta x(t-\tau ).$$

(177)

When the solution to Equation (177) is bounded, the solution to Equation (166) is Lyapunov stable. When the solution to Equation (177) converges to zero for $t\to \infty$, a fixed point is asymptotically stable.

The trivial fixed point, $x_1^{*}=0$, is stable when

$$\mathrm{sgn}\gamma =-1,\mathrm{sgn}B=\pm 1,-\infty <b<\infty ,\tau \ge 0,$$

(178)

but it is always unstable for $\mathrm{sgn}\gamma =1$,

The motion in the vicinity of nontrivial fixed points is described by Equation (177), which, using relation (169), reduces to the equation,

$$\frac{d}{dt}\delta x\left(t\right)=-\left(\mathrm{sgn}\gamma \right)\delta x\left(t\right)+b\left(\mathrm{sgn}\gamma \right)x^{*}\delta x(t-\tau ).$$

(179)

Looking for the solution in the form,

$$\delta x\left(t\right)\propto e^{-\lambda t},$$

(180)

one gets the equation for the Lyapunov exponent $\lambda$ depending on the studied fixed point.

3.3. Punctuated Evolution

In the biological evolution theory there exists a hypothesis, called punctuated equilibrium, suggesting that the evolutional changes of biological species are marked by episodes of rapid speciation between long periods of little or no change. This type of evolution, that occurs rapidly, being separated by periods of stasis, or equilibrium, is called punctuated equilibrium [146,147,148]. If biological changes can be described by a quantitative characteristic, then the corresponding graph has the shape of a ladder. A mathematical model of such a punctuated development can be represented by a delayed equation [145]. This model might characterise not only the evolution of biological species, but also the evolution of firms and other organisations. In order to emphasise that this type of development is rather general, and can occur not only for biological species, it is called punctuated evolution.

Let us consider the most realistic case of population characterised by gain and competition [145], where

$$\mathrm{sgn}\gamma =\mathrm{sgn}B=1.$$

(181)

Then, Equation (166) reads:

$$\frac{dx}{dt}=x-x^{2}exp\{-bx(t-\tau )\}.$$

(182)

In that case, the trivial fixed point, $x_1^{*}=0$, is never stable. The nontrivial fixed point, $x_2^{*}$, is stable when either

$$0<x_2^{*}\le \frac{1}{e}(b\le -e,\tau <{\tau }_2),$$

(183)

where

$${\tau }_2=\frac{1}{\sqrt{{\left(b{x}_2^{*}\right)}^2-1}}arccos\left(\frac{1}{b{x}_2^{*}}\right),$$

(184)

or in the range,

$$\frac{1}{e}<x_2^{*}<e\left(-e<b\le \frac{1}{e},\tau >0\right).$$

(185)

The fixed point, $x_3^{*}>e$, is always unstable under gain and competition. So that the sole stable fixed point is $x_2^{*}\in (0,e)$.

When the production parameter, b, is semi-negative, which implies that the population does not produce its carrying capacity but rather destroys it or, in the best case, retains the given capacity value, then the population remains bounded. More precisely, the following theorem takes place [145].

Theorem on population boundedness. The solution $x\left(t\right)$ to the evolution Equation (182), for $b\le 0$, any finite $\tau \ge 0$, and any initial conditions $x_0\ge 0$, is bounded for all times $t\ge 0$, and, for $b<0$, there exists a time $t_0=t_0(x_0,\tau )$ such that

$$0\le x\left(t\right)\le 1(t>t_0).$$

The proof is given in Ref. [145].

It is important to note that, for some values of the production parameter b, the basin of attraction of $x_2^{*}$ is not the whole positive semi-line $x_0\ge 0$, but a limited interval. This happens for $b\in (0,1/e)$, if the basin of attraction is given by the inequalities

$$0<x_0<x_3^{*}\left(0<b\le \frac{1}{e}\right).$$

(186)

For the values of b satisfying condition (186), but with $x_0>x_3^{*}$, the solution x tends to infinity, as $t\to \infty$.

Overall, there exist the following regimes of population dynamics:

(i) Punctuated unbounded growth. This growth happens when b is outside of the stability region of $x_2^{*}$, so that

$$b>\frac{1}{e},\tau >0,x_0>0.$$

(187)

An analogous unbounded punctuated growth happens when b is inside the stability region, but $x_0$ is outside of the attraction basin of $x_2^{*}$, which occurs for

$$0<b<\frac{1}{e},\tau >0,x_0>x_3^{*}.$$

(188)

The unbounded punctuated behaviour happens when either the production parameter is sufficiently large or the initial level of creative activity is enough high.

(ii) Punctuated convergence to a bounded state. Punctuated convergence to a finite stationary state $x_2^{*}$ happens for positive production parameters, when

$$0<b<\frac{1}{e},\tau >0,x_0<x_3^{*}.$$

(189)

The solutions tend to the stationary state $x_2^{*}$ by punctuated steps from below, if $x_0<x_2^{*}$, and from above, if $x_0>x_2^{*}$.

(iii) Oscillatory convergence to a bounded state. When the production parameter is negative, the approach to the stationary state becomes oscillatory. For

$$-e<b<0,\tau \ge 0,x_0>0,$$

(190)

there appear sharp reversals after almost horizontal plateaus. For

$$b<-e,\tau \le {\tau }_2,x_0>0,$$

(191)

the population dynamics becomes strongly oscillating when approaching to the focus $x_2^{*}$.

(iv) Everlasting oscillations. When the destructive action is rather strong and the time delay is long, so that

$$b<-e,\tau \ge {\tau }_2,x_0>0,$$

(192)

then there develops the regime of everlasting oscillations.

The discussed regimes of dynamics under gain and competition are described in detail in Ref. [145], as well as other cases corresponding to gain and cooperation, loss and competition, and loss and cooperation. The punctuated evolution is typical only for the regime of gain and competition. This most interesting regime is rather widespread for many population societies. Thus it occurs in the evolution of biological species, growth of world population and urban population, development of science and technology, development of art and culture, energy production and consumption, growth of social organisations, growth of alive organisms, the improvement of individual abilities, proliferation of cells and bacteria, as well as in the decay of biological organisms and societies [123,141,145,146].

3.4. Symbiosis of Species

The term “symbiosis” characterises close and sufficiently long-time interactions between different biological species. In biological and ecological societies, symbiotic relationships are widespread. For example, the symbiosis between plant roots and fungi is a typical feature of numerous ecosystems. Co-existence of coral reefs and fishes is another well known example. Numerous other examples can be found elsewhere [149,150,151,152,153]. In human societies, the examples of symbiosis are also widespread. Symbiotic relations are common for firms and banks, people and government, culture and language, economic and intellectual levels of society, basic science and technological applications.

3.4.1. Interaction through Carrying Capacities

The main idea in the new model of symbiosis is the observation that in symbiotic relations it is not the species themselves that interact directly with each other, but that symbiotic species influence the carrying capacities of each other. This implies that in the evolution equations of species the carrying capacities, $K_j$, are functions of the populations $N_1$, $N_2$, etc. It is convenient to normalise the populations $N_i$ with respect to different normalisation constants, $M_i$, introducing the fractions

$$x_i\equiv \frac{N_i}{M_i}(M_i=\mathrm{const}).$$

(193)

Then the reduced form of the equations, without external forces, can be written as

$$\frac{d{x}_i}{dt}=\left({\gamma }_i-\sum _j\frac{B_{ij}{M}_j}{K_j}{x}_j\right)x_i.$$

(194)

The convenience of using different normalisations is easily understood if one remembers that different species can exhibit rather different numbers of their members. For instance, the most important to humans symbiosis is that one between the human body and the numerous organisms of the microbiome, where there are more than 2000 bacterial species and about ${10}^{14}$ microorganisms. The latter number is about ten times the number of cells of the human body [154]. Essentially different numbers of different populations can require different normalisation constants.

The coexisting species are characterised by direct interactions, $B_{ij}$, and the interactions through their carrying capacities. In the general case, one can take into account the whole matrix of the interaction elements $B_{ij}$. When the direct interactions between the members of the same kind of species are much larger than between the members of different species, then the matrix $B_{ij}$ is approximately diagonal,

$$B_{ij}\approx B_{ii}{\delta }_{ij}.$$

(195)

Therefore, Equation (194) reduces to

$$\frac{d{x}_i}{dt}={\gamma }_i\left(1-\frac{B_{ii}{M}_i}{{\gamma }_i{K}_i}{x}_i\right)x_i.$$

(196)

It is possible to introduce dimensionless carrying capacities,

$$q_i\equiv \frac{K_i}{K_i^0},$$

(197)

that are functions $q_i=q_i(x_1,x_2,\dots )$ such that $q(0,0,\dots )=1$, which is possible by defining the appropriate normalisation constants $K_i^0$. Generally, the variables entering $q_i$ should be delayed in time, such that $x_i=x_i(t-{\tau }_i)$.

Let us consider the case when the parameters ${\gamma }_i$ and $B_{ii}$ are positive. And let us choose the normalisation factors, $M_i$, as

$$M_i\equiv K_i^0\frac{{\gamma }_i}{B_{ii}}.$$

(198)

Then,

$$\frac{d{x}_i}{dt}={\gamma }_i\left(1-\frac{x_i}{q_i}\right)x_i.$$

(199)

Considering the case of two species, $x_1\equiv x$ and $x_2\equiv y$, Equation (199) reduces to the system of two equations

$$\frac{dx}{dt}={\gamma }_1\left(1-\frac{x}{q_1}\right)x,\frac{dy}{dt}={\gamma }_2\left(1-\frac{y}{q_2}\right)x.$$

(200)

Let us measure time in units of $1/{\gamma }_2$ and define the ratio

$$\alpha \equiv \frac{{\gamma }_1}{{\gamma }_2}.$$

(201)

Thus, one comes to the equations,

$$\frac{dx}{dt}=\alpha \left(x-\frac{x^2}{q_1}\right),\frac{dy}{dt}=y-\frac{y^2}{q_2}.$$

(202)

Now, one needs to define the carrying capacities $q_i$ that, generally, are the functions of the delayed populations,

$$q_i=q_i(x_1(t-{\tau }_1),x_2(t-{\tau }_2),\dots ).$$

(203)

For the case of two populations, $q_i=q_i(x,y)$. If the function $q_i(x,y)$ to be expanded in a series in powers of populations, one gets it in the form as

$$q_i\simeq 1+c_{1}x+c_{2}y+c_{11}{x}^2+c_{22}{y}^2+c_{12}xy+\dots .$$

(204)

If the mutual influence of populations on the carrying capacities of each other is weak, then $q_i$ can be expanded over populations, with limiting ourselves by the lowest terms of the expansion [123]. In the case of strong influence of symbiotic populations on each other, limiting ourselves by several first terms gives an expression that can become zero, thus producing spurious divergences in the terms containing $1/q_i$. To include in the consideration a strong mutual influence, the expansions can be summed by means of exponential summation [142,143,144], thus avoiding spurious zeroes in the effective carrying capacity [155,156]. In that way, one can derive the form,

$$q_i=exp\left\{\sum _j{b}_{ij}{x}_j+\sum _{jk}{c}_{ijk}{x}_j{x}_k\right\}.$$

(205)

The first term in the exponential describes the influence on the carrying capacity of separate populations not correlated with each other. The second term in the exponential characterises the impact of the populations on the carrying capacity, when the symbiotic populations correlate with each other.

From the general expression (205), it is possible to set off two typical cases. One case is when the main terms in the exponential are those corresponding to the action of the symbiotic species without their mutual correlations, that is when the carrying capacities, for the case of two species, have the form

$$q_1=e^{by},q_2=e^{gx}\left(\mathrm{uncorrelated}\mathrm{symbiosis}\right),$$

(206)

which can be called uncorrelated symbiosis, since the species, in the process of their action on the carrying capacities, do not correlate with each other. And the other scenario describes correlated symbiosis, when the species influence the carrying capacities being mutually correlated, which for two species reads:

$$q_1=e^{bxy},q_2=e^{gxy}\left(\mathrm{correlated}\mathrm{symbiosis}\right).$$

(207)

The self-action of the species on their own carrying capacities is neglected here assuming that the influence of the symbiotic species is more important. It is also possible to study a mixed case when $q_1$ corresponds to a correlated symbiosis, while $q_2$, to uncorrelated symbiosis.

Representing symbiosis through the mutual action of the symbiotic species on the carrying capacities of each other allows for the description of all types of symbiosis, which is straightforwardly connected with the signs of the symbiotic parameters b and g. Below, the classification of symbiotic types is given. In order to be precise and not to disturb the meaning, the formulation of the definitions below are given closely following Refs. [155,156].

(i) Mutualism is a relationship between different species when both of them receive mutual benefit, which corresponds to the symbiotic parameters,

$$b>0,g>0.$$

(208)

(ii) Commensalism is a relationship, when one of the species benefits from the coexistence with the other species, while the other one is neutral, getting neither profit nor harm, which is defined by one of the conditions,

$$b>0,g=0,$$

(209)

$$b=0,g>0.$$

(210)

(iii) Parasitism is a relation, when at least one of the coexisting species is harmful to the other one, which is characterised by one of the conditions,

$$b>0,g<0,$$

(211)

$$b=0,g<0,$$

(212)

$$b<0,g<0,$$

(213)

$$b<0,g=0,$$

(214)

$$b<0,g>0.$$

(215)

As is seen, there exist various types of symbiosis described by the systems of differential equations that can be solved numerically [157]. Since symbiosis is extremely widespread in nature, there are numerous particular examples of coexisting species among biological and ecological societies, including bacteria and viruses [149,150,151,152,153,158].

3.4.2. Uncorrelated Symbiosis

Dynamics of symbiotic species, under the uncorrelated symbiosis and the assumption of approximately equal rates ${\gamma }_1$ and ${\gamma }_2$, are described by the system of equations,

$$\frac{dx}{dt}=x-x^2{e}^{-by},\frac{dy}{dt}=y-y^2{e}^{-gx},$$

(216)

with the symbiotic parameters $b\in (\infty ,\infty )$ and $g\in (-\infty ,\infty )$. For simplicity, the delay times of the populations entering the carrying capacities are assumed considerably small, so that can be neglected. The initial conditions are $x_0=x\left(0\right)$ and $y_0=y\left(0\right)$.

For any values of the parameters, there always exist three trivial fixed points, $\{x^{*}=0,y^{*}=0\}$, $\{x^{*}=1,y^{*}=0\}$, and $\{x^{*}=0,y^{*}=1\}$, which are unstable for all g and b. Nontrivial fixed points are given by the equations,

$$x^{*}=e^{b{y}^{*}},y^{*}=e^{g{x}^{*}},$$

(217)

that can be rewritten as

$$x^{*}=exp\left(b{e}^{g{x}^{*}}\right),y^{*}=exp\left(g{e}^{b{y}^{*}}\right).$$

(218)

The Lyapunov exponents are defined by the expressions

$${\lambda }_1=-1+\sqrt{bg{x}^{*}{y}^{*}},{\lambda }_2=-1-\sqrt{bg{x}^{*}{y}^{*}}.$$

(219)

Under mutualism, the analysis shows that there can exist the following regimes depending on the parameters $b>0$ and $g>0$:

(i)

unbounded growth of populations with time,

(ii)

convergence to a stationary state.

In the case of parasitism, the picture depends on whether a single species is parasitic or both species are parasites. When a single species is parasitic, that is when either $b<0$, while $g>0$ or $b>0$, while $g<0$, then only one regime exists, when the populations tend to their stationary states.

When both species are parasitic, so that $b<0$ and $g<0$ then, depending on the symbiotic parameters, there can exist two regimes:

(i)

convergence to single stationary state,

(ii)

bistability with two stationary states;

see Ref. [156] for details.

3.4.3. Correlated Symbiosis

Correlated symbiosis of two species is characterised by the system of equations,

$$\frac{dx}{dt}=x-x^2{e}^{-bxy},\frac{dy}{dt}=y-y^2{e}^{-gxy},$$

(220)

where it is again, as above, assumed that ${\gamma }_1$ is close to ${\gamma }_2$. Similarly to the case of uncorrelated symbiosis, there always exist three trivial fixed points, $\{x^{*}=0,y^{*}=0\}$, $\{x^{*}=1,y^{*}=0\}$, and $\{x^{*}=0,y^{*}=1\}$, which are unstable for any g and b. Nontrivial fixed points are the solutions to the equations,

$$x^{*}=e^{b{x}^{*}{y}^{*}},y^{*}=e^{g{x}^{*}{y}^{*}},$$

(221)

that can be rewritten as

$$x^{*}=exp\left\{b{\left(x^{*}\right)}^{1+g/b}\right\},y^{*}=exp\left\{g{\left(y^{*}\right)}^{1+b/g}\right\}.$$

(222)

The Lyapunov exponents (36) are defined by the relations,

$${\lambda }_1=-1,{\lambda }_2=-1+(b+g)x^{*}{y}^{*}.$$

(223)

In the case of mutualism, where $b>0$ and $g>0$, there are the following types of behaviour:

(i)

unbounded growth of both populations,

(ii)

convergence of populations to stationary states.

In the case of parasitism, when either one of the symbiotic parameters is negative, while the other one is positive, or both parameters are negative, the following scenarios are possible:

(i)

convergence to stationary states,

(ii)

unbounded growth of parasitic population and dying out of host population;

see Ref. [156] for detailed investigation.

3.4.4. Mixed Symbiosis

In the examples describing symbiosis, as given in Section 3.4.1, Section 3.4.2 and Section 3.4.3, it is possible to notice the symmetry corresponding to the interchange between the populations x and y. It also may happen that two symbiotic species exhibit non-symmetric relations that can be characterised by the case of mixed symbiosis, when one of the species displays correlated symbiosis, while the other species, uncorrelated symbiosis, This case is described by the equations,

$$\frac{dx}{dt}=x-x^2{e}^{-bxy},\frac{dy}{dt}=y-y^2{e}^{-gx}.$$

(224)

Similarly to the cases, dicussed above, there always exist three trivial fixed points, $\{x^{*}=0,y^{*}=0\}$, $\{x^{*}=1,y^{*}=0\}$, and $\{x^{*}=0,y^{*}=1\}$, which are unstable for all symbiotic parameters b and g. Nontrivial fixed points are given by the equations,

$$x^{*}=e^{b{x}^{*}{y}^{*}},y^{*}=e^{g{x}^{*}}$$

(225)

that can be represented as

$$x^{*}=exp\left\{b{x}^{*}{e}^{g{x}^{*}}\right\},y^{*}=exp\left\{g{\left(y^{*}\right)}^{b{y}^{*}/g}\right\}.$$

(226)

The Lyapunov exponents (36) are:

$$\begin{gathered} {\lambda }_1=\frac{1}{2}\left[b{x}^{*}{y}^{*}-2+x^{*}\sqrt{b{y}^{*}(4g+b{y}^{*})}\right], \\ {\lambda }_2=\frac{1}{2}\left[b{x}^{*}{y}^{*}-2-x^{*}\sqrt{b{y}^{*}(4g+b{y}^{*})}\right]. \end{gathered}$$

(227)

The dynamics of the symbiotic populations strongly depends on whether the influence of the species y on the carrying capacity of species x is mutualistic or parasitic, which is described by the sign of the symbiotic parameter b. If the latter is negative $(b<0)$, there exists the sole regime, when both populations tend to their stationary limits. If $b>0$ and $g>0$, two regimes can happen:

(i)

unbounded growth of both populations,

(ii)

convergence to stationary states.

When $b>0$, but the first species is parasitic, with $g<0$, then there can exist two possibilities:

(i)

convergence to stationary states,

(ii)

everlasting oscillations;

see Ref. [156] for details. Applications to financial markets, considering symbiosis between asset prices and bonds, resulting in periodically growing and collapsing bubbles, are analysed in Ref. [159].

3.5. Role of Growth Rates

When the growth rates ${\gamma }_1$ and ${\gamma }_2$ are essentially different, it is necessary to consider the system of equations,

$$\frac{dx}{dt}=\alpha \left(x-\frac{x^2}{q_1}\right),\frac{dy}{dt}=y-\frac{y^2}{q_2},$$

(228)

in which $\alpha \equiv {\gamma }_1/{\gamma }_2$ and time is measured in units of ${\gamma }_2^{-1}$. Without the loss of generality, it is possible to take $\alpha >1$. In the case of uncorrelated symbiosis, one has:

$$\frac{dx}{dt}=\alpha \left(x-x^2{e}^{-by}\right),\frac{dy}{dt}=y-y^2{e}^{-gx}.$$

(229)

There are three trivial fixed points: the unstable node $\{0,0\}$, with the Lyapunov exponents, ${\lambda }_1=1$ and ${\lambda }_2=\alpha$; a saddle $\{1,0\}$, with the Lyapunov exponents, ${\lambda }_1=1$ and ${\lambda }_2=-\alpha$; and the saddle $\{0,1\}$, with the Lyapunov exponents ${\lambda }_1=-1$ and ${\lambda }_2=\alpha$. The nontrivial fixed points do not depend on the value of $\alpha$ and are defined as in Section 3.4.2. The related Lyapunov exponents are:

$$\begin{gathered} {\lambda }_1=-\frac{1}{2}(1+\alpha )+\sqrt{{(1-\alpha )}^2+4\alpha bg{x}^{*}{y}^{*}}, \\ {\lambda }_2=-\frac{1}{2}(1+\alpha )-\sqrt{{(1-\alpha )}^2+4\alpha bg{x}^{*}{y}^{*}}. \end{gathered}$$

(230)

Although the fixed points do not depend on the rate $\alpha$, but the Lyapunov exponents do depend on it, as well as the related basins of attraction also depend on $\alpha$, which hence influences the stability of fixed points [160]. Therefore by varying $\alpha$ it is admissible to obtain different regimes of motion, hence to realise dynamic phase transitions by the sole variation of the growth rate, with keeping other parameters unchanged.

A dynamic transition implies a qualitative change of dynamical behaviour, when system parameters reach a bifurcation point. Then the type of fixed points changes. However, there can happen a nonstandard dynamic transition, when a qualitative change of dynamical behaviour occurs due to the variation of the growth rate. In that case, the fixed points do not change, remaining the same, while a sharp change in dynamical behaviour happens because the growth rate shifts the boundary of the basins of attraction, so that the initial point of a trajectory, which was inside the attraction basin, moves outside of the basin [161].

3.6. Self-Organised Society

The evolution of biological societies, including human societies, is an especially important old problem studied in voluminous literature beginning with Darwin [162,163]. Societies are typically structured into groups representing particular features or traits. As an example of group representatives, it is possible to mention collaborators and defectors. The collaborators cooperate with each other for the benefits of the whole society, while defectors, on the contrary, exploit it [101,102,103].

The evolution of groups is commonly studied on the basis of the replicator equation discussed in the above sections. As has been explained at the beginning of this Section, if the society consists of two groups only, cooperators and defectors, the latter always outperform the former, so that the sole evolutionary stable state is the state where there are no cooperators, but there exist solely defectors. It is understandable that for a closed self-organised society, where there are no unlimited resources supplied from somewhere outside, this conclusion is absurd, since defectors produce nothing and, being left alone, cannot survive.

The aim here is to consider a self-organised society, whose means of survival are produced inside the society itself. A closed self-organised society cannot exist being composed solely of defectors because they will have no means for survival. Sometimes, to correct this unrealistic conclusion, one introduces punishers. However, these also require means for their existence and cannot survive if nothing is produced.

3.6.1. Trait Groups

In this Section, a self-organised structured society is described composed of the trait groups representing four types of typical agents: cooperators, defectors, regulators, and outsiders. The main novelty of the approach is that, instead of a single fitness or utility for each group, relative utilities for each group with respect to the society as a whole, and mutual utilities with respect to each other are introduced [164].

Let us consider a society composed of several groups, whose fractions are defined as the ratios,

$$x_i\equiv \frac{N_i}{N_0},$$

(231)

with respect to a normalising number $N_0$ that can be chosen to define the initial total number of the society members,

$$N_0=\sum _i{N}_i\left(0\right)(t=0).$$

(232)

Then, at the initial time there exists the normalisation,

$$\sum _i{x}_i\left(0\right)=1.$$

(233)

However, at the later times such a normalisation, generally, is not valid. In general, the number of the society members can vary with time because of births, deaths, and the influx of outsiders.

Biological, including human, societies have much in common with the structures of the organisms because of which the society groups can be straightforwardly compared with the parts of biological organisms. The classification below follows Refs. [165,166].

• Cooperators $\left(x_1\right)$, who contribute to the whole society. In a human society, the cooperators form the working force producing the gross domestic product. In a biological organism, cooperators can be associated with healthy cells.
• Defectors $\left(x_2\right)$, who do not contribute to the society and can exist only owing to the work of cooperators. In a social system, the groups that benefit from the society support without contributing are prisoners, pensioners, and unemployed people. In a biological organism, defectors can be represented by ill cells.
• Regulators $\left(x_3\right)$, who maintain order in the society and punish defectors and harmful outsiders. In a human society, this role is played by the police, the army, and the order enforcing bureaucracy. To support the existence of regulators, the society has to pay the necessary costs. In a biological organism, regulators can correspond to the cells of the immune system.
• Outsiders $\left(x_4\right)$, who also exploit the society, but, contrary to defectors, the difference is that they enter the society from outside. The harmful outsiders could be interpreted as terrorists or as foreign invading armies. For biological organisms, outsiders could be pathogens or viruses infecting the organism.

As the evolution equations for fractions (231), it is possible to accept the rate equations,

$$\frac{d{x}_i}{dt}=\left({\gamma }_i+\sum _j{a}_{ij}{x}_j\right)x_i+{\phi }_i.$$

(234)

These equations are formally equivalent to the rate Equation (52) or to the Lotka–Volterra Equation (93). However, the meaning of the parameters entering Equation (234) is different. In the Lotka–Volterra equations, ${\gamma }_i$ are the birth-death rates and $a_{ij}$ are intensities of interactions.

In the case here considred, each ${\gamma }_i$ plays the role of utility rate, or the rate of production (or consumption), or the production of resources with respect to the whole society, and the parameters $a_{ij}$ are interpreted as the utility production (or consumption) by the jth group with respect to the ith group. Diagonal elements $a_{ii}$ correspond to the competition between the members of the same group (hence $a_{ii}<0$). The signs of non-diagonal $a_{ij}$ are defined depending on the usefulness of the group j to group i, so that if j is useful to i, then $a_{ij}$ is positive, while, when j is harmful for i, then $a_{ij}$ is negative. External influx is assumed to exist only for outsiders. The definitions below are based on Refs. [165,166].

Cooperators are useful for the society. Defectors are not useful for cooperators. Regulators require the support of cooperators which is a cost to the latter. Harmful insiders also are not useful for cooperators. Summarising this, one has:

$${\gamma }_1>0,{\phi }_1=0,a_{11}<0,a_{12}<0,a_{13}<0,a_{14}<0.$$

(235)

Defectors are not useful for the society. Cooperators are necessary for defectors who live at the expense of the former. Regulators suppress and punish defectors. Invaders are not useful for defectors. Thus:

$${\gamma }_2<0,{\phi }_2=0,a_{21}>0,a_{22}<0,a_{23}<0,a_{24}<0.$$

(236)

Regulators do not produce resources. Society needs to support regulators. Cooperators are necessary for regulators. The role of regulators is to maintain order and to punish defectors, whose presence justifies the existence of regulators. Similarly, regulators suppress harmful invaders, which justifies the existence of regulators. Therefore:

$${\gamma }_3<0,{\phi }_3=0,a_{31}>0,a_{32}>0,a_{33}<0,a_{34}>0.$$

(237)

Outsiders are not useful for society. However, cooperators are necessary for outsiders. Outsiders exploit defectors by taking a part of their share, and benefit from their presence. Regulators, suppressing outsiders, are not useful to them. Hence:

$${\gamma }_4<0,{\phi }_4\ge 0,a_{41}>0,a_{42}>0,a_{43}<0,a_{44}<0.$$

(238)

Taking into account the signs of the parameters, Equation (234) is written as

$$\begin{gathered} \frac{d{x}_1}{dt}=({\gamma }_1-|a_{11}|x_1-|a_{12}|x_2-|a_{13}|x_3-|a_{14}|x_4)x_1, \\ \frac{d{x}_2}{dt}=(-|{\gamma }_2|+|a_{21}|x_1-|a_{22}|x_2-|a_{23}|x_3-|a_{24}|x_4)x_2, \\ \frac{d{x}_3}{dt}=(-|{\gamma }_3|+|a_{31}|x_1+|a_{32}|x_2-|a_{33}|x_3+|a_{34}|x_4)x_3, \\ \frac{d{x}_4}{dt}=(-|{\gamma }_4|+|a_{41}|x_1+|a_{42}|x_2-|a_{43}|x_3-|a_{44}|x_4)x_4+{\phi }_4. \end{gathered}$$

(239)

In general, the quantities $a_{ij}$ could be functions of populations, as in the considered above cases of retarded carrying capacity, resulting in punctuated evolution, and symbiotic relations through functional carrying capacities. However, not to complicate the consideration, $a_{ij}$ to be treated as parameters. To reduce the number of parameters, it is possible to employ relations existing between them, keeping in mind that, by accepted interpretation, $a_{ij}$ is the utility rate for a group i of a group j. For brief, below $a_{ij}$ are called utilities.

The cooperators are the sole group producing resources for the whole society. These resources are denoted through ${\gamma }_1$. The cooperators compete for these resources, which is described by the value $|a_{11}|$. This utility is actually all that is produced by the cooperators, which implies the equality,

$$|a_{11}|={\gamma }_1.$$

(240)

The diagonal elements $|a_{ii}|$ describe the strength of competition among the members of a group i for the available for them resources ${\gamma }_i$ playing the role of a carrying capacity. The standard dependence for the term characterising competition on the carrying capacity is the inverse dependence,

$$|a_{ii}|=\frac{C}{|{\gamma }_i|},$$

(241)

where C is a constant. From relation (240) it follows that $C={\gamma }_1^2$. Therefore

$$|a_{ii}|=\frac{{\gamma }_1^2}{|{\gamma }_i|}.$$

(242)

The non-diagonal elements $|a_{ij}|$ play the role of mutual utility for the groups i and j. The widely used expression for the mutual utility is the Bernoulli-Nash mutual utility [167,168] that can be written as

$$|a_{ij}|=\sqrt{{\gamma }_i{\gamma }_j}(i\ne j).$$

(243)

Additionally, there should exist the reciprocity relation,

$$a_{ij}=-a_{ji}(i\ne j)$$

(244)

signifying that, if a group i receives some utility from a group j, then the group j loses the same amount of utility.

Let us introduce the dimensionless parameters, describing the fractions of the wealth consumed by the related groups with respect to the total amount of resources produced by the society cooperators,

$$a\equiv \frac{|{\gamma }_2|}{{\gamma }_1},b\equiv \frac{|{\gamma }_3|}{{\gamma }_1},c\equiv \frac{|{\gamma }_4|}{{\gamma }_1},$$

(245)

and the dimensionless influx of outsiders

$$\phi \equiv \frac{{\phi }_4}{{\gamma }_1}.$$

(246)

If time is measured in units of ${\gamma }_1^{-1}$, then Equation (239) acquires the form of the system of equations

$$\frac{d{x}_i}{dt}=f_i(i=1,2,3,4),$$

(247)

with the right-hand sides,

$$\begin{gathered} f_1=\left(1-x_1-\sqrt{a}{x}_2-\sqrt{b}{x}_3-\sqrt{c}{x}_4\right)x_1, \\ {f}_2=\left(-a+\sqrt{a}{x}_1-\frac{1}{2}{x}_2-\sqrt{ab}{x}_3-\sqrt{ac}{x}_4\right)x_2, \\ {f}_3=\left(-b+\sqrt{b}{x}_1+\sqrt{ab}{x}_2-\frac{1}{b}{x}_3+\sqrt{bc}{x}_4\right)x_3, \\ {f}_4=\left(-c+\sqrt{c}{x}_1+\sqrt{ac}{x}_2-\sqrt{bc}{x}_3-\frac{1}{c}{x}_4\right)x_4. \end{gathered}$$

(248)

One more reasonable restriction that should be accepted is the conservation law telling that one cannot consume more than it is produced. That is, what is consumed by defectors, regulators, and outsiders, cannot be larger than what is produced by cooperators,

$$|{\gamma }_2+{\gamma }_3+{\gamma }_4|\le {\gamma }_1.$$

(249)

In dimensionless units this is equivalent to the inequality,

$$a+b+c\le 1.$$

(250)

If the condition (250) becomes invalid, this means that defectors, regulators, and outsiders consume more than cooperators produce. This could be realised only under a supply from outside or by the destruction of the given carrying capacity, which implies instability.

3.6.2. Coexistence of Cooperators and Defectors

It is instructive to consider the coexistence of two typical groups, cooperators and defectors. Recall that in the case of replicator equation, the stable state corresponds to the absence of cooperators and all society being composed of defectors. As has been discussed, this absolutely unrealistic conclusion stems from the fact that the replicator equation does not describe a self-organised society. Let us see what happens in the case of a self-organised society.

Considering cooperators and defectors, one has the evolution equations,

$$\frac{d{x}_1}{dt}=\left(1-x_1-\sqrt{a}{x}_2\right)x_1,\frac{d{x}_2}{dt}=\left(-a+\sqrt{a}{x}_1-\frac{1}{a}{x}_2\right)x_2.$$

(251)

Looking for fixed points satisfying conditions (250), one finds the only solution,

$$x_1^{*}=\frac{1+a^{5/2}}{1+a^2},x_2^{*}=\frac{1-\sqrt{a}}{1+a^2}{a}^{3/2},$$

(252)

that is stable under the condition,

$$0\le a\le 1.$$

(253)

From Equation (252), one finds that the minimal fraction of cooperators occurs when

$$\underset{a}{min}{x}_1^{*}=0.940(a=0.565),$$

(254)

while the maximal fraction of defectors is

$$\underset{a}{max}{x}_2^{*}=0.083(a=0.480).$$

(255)

Thus in a stable self-organised society, the amount of defectors cannot be larger than around $10\%$. If this amount is essentially larger, the society is not stable.

3.6.3. Three Coexisting Groups

The standard picture is the coexistence of three groups, cooperators, defectors, and regulators, while outsiders are not numerous, hence can be neglected. Then, the equations are:

$$\begin{gathered} \frac{d{x}_1}{dt}=\left(1-x_1-\sqrt{a}{x}_2-\sqrt{b}{x}_3\right)x_1, \\ \frac{d{x}_2}{dt}=\left(-a+\sqrt{a}{x}_1-\frac{1}{a}{x}_2-\sqrt{ab}{x}_3\right)x_2, \\ \frac{d{x}_3}{dt}=\left(-b+\sqrt{b}{x}_1+\sqrt{ab}{x}_2-\frac{1}{b}{x}_3\right)x_3. \end{gathered}$$

(256)

There is the only fixed point satisfying the conservation law (250),

$$\begin{gathered} x_1^{*}=\frac{1+a^2{b}^2+a^{5/2}(1+b^2)+b^{5/2}(1-a^2)}{(1+a^2)(1+b^2)}, \\ {x}_2^{*}=\frac{1-\sqrt{a}(1+b^2)-b^2(1-2\sqrt{b})}{(1+a^2)(1+b^2)}{a}^{3/2},x_3^{*}=\frac{1-\sqrt{b}}{1+b^2}{b}^{3/2}, \end{gathered}$$

(257)

which is stable when

$$0\le a<{\left(\frac{2b^2\sqrt{b}-b^2+1}{1+b^2}\right)}^2,0\le b\le 1.$$

(258)

This shows that the minimal fraction of cooperators is

$$\underset{a,b}{min}{x}_1^{*}\approx 0.91(a\approx 0.45,b\approx 0.55),$$

(259)

the maximal fraction of defectors is

$$\underset{a,b}{max}{x}_2^{*}\approx 0.08(a\approx 0.48,b=0),$$

(260)

and the maximal fraction of regulators is

$$\underset{a,b}{max}{x}_3^{*}\approx 0.08(0\le a\le 0.52,b=0.48).$$

(261)

Again, as in Section 3.6.2, one can see that in a stable society the fractions of defectors and regulators should not exceed about $10\%$ each.

More discussions and applications of the theory for describing ant and bee colonies, and concrete countries is given in Ref. [164], where the influence of noise is also studied.

4. Models of Financial Markets

Financial markets can be treated as complex social systems, because of which their behaviour could allow for mathematical description. It is not our aim to plunge deeply into the ocean of economic theories, but our aim is to present some basic ideas the models of markets are based on and to illustrate them by simple enough examples.

4.1. Efficient Market Model

The old standing hypothesis that the functioning of financial markets is based on, is the efficient market hypothesis. It is generally accepted [169,170] that the idea of efficient market has been anticipated by Bachelier [171], who compared the motion of market prices with the Brownian motion. According to Bachelier, “past, present and even discounted future events are reflected in market price, but often show no apparent relation to price changes”. Fama [172,173] distinguished three forms of efficient market hypothesis, weak form, semi-strong form, and strong form. Weak form assumes that the stock prices indicate the present public market information, the past performance does not play any role, although the prices may not reflect new information that is not yet publicly available. Semi-strong form states that the stock prices reflect both the market and non-market public information and that prices adjust quickly to any new public information that becomes available. Strong form insists that prices reflect the entirety of both public and private information, historical and new, as well as insider information. Summarising, the strong form of the efficient market hypothesis assumes that:

• all prices on traded assets already reflect all past available information;
• all prices instantly adjust to any newly appearing information;
• all prices instantly reflect even hidden information.

An efficient market is assumed to be absolutely equilibrium, random, and stable. It unpredictably fluctuates, so that nobody can consistently achieve returns in excess of average market returns. The motion of stock prices is compared with random walk typical of Brownian motion.

Typical agents acting in a market are supposed to be rational. Any one particular agent can be wrong and make mistakes, but the market as a whole is always right. Agents errors are random, so that they are averaged out. A typical, or representative agent is the so-called Homo Economicus who:

• possesses all existing information necessary for the correct price evaluation;
• can immediately process all existing information;
• makes objective unbiased conclusions based on the maximisation of expected utility.

Let $p=p\left(t\right)$ represent asset price, stock price, or market index. Dynamics of p, for an efficient market, is assumed to follow the geometric Brownian motion that satisfies the stochastic differential equation.

$$dp\left(t\right)=[ydt+\sigma dW(t\left)\right]p\left(t\right),$$

(262)

in which y is a drift rate of the price, $\sigma$ the diffusion coefficient, or market volatility of returns, and W is a random Wiener process. Note that the above equation can be interpreted as a rate equation,

$$\frac{dp\left(t\right)}{dt}=R\left(t\right)p\left(t\right),$$

(263)

with a randomly fluctuating rate,

$$R\left(t\right)=y+\sigma \xi \left(t\right),$$

(264)

where $\xi \left(t\right)=dW\left(t\right)/dt$ is white noise. This model describes the market motion as on average exponential, with superimposed random fluctuations.

Without the noise, the solution would be

$$p\left(t\right)=p\left(0\right)e^{yt}.$$

(265)

While in the presence of noise, the numerical solution is given by the iterative scheme,

$$p\left(t_{n+1}\right)=[1+y\Delta t+\sigma R\left(t_n\right)\sqrt{\Delta t}]p\left(t_n\right).$$

(266)

4.2. Diffusion Price Model

In the diffusion price model, the drift term y in Equation (262) is treated not as a constant, but as a random variable, so that the market is described by the system of equations,

$$dp=(ydt+{\sigma }_{1}d{W}_1)p,dy=-\gamma ydt+{\sigma }_{2}d{W}_2,$$

(267)

where $\gamma >0$ is called “market friction”. The initial conditions are $p\left(0\right)$ and $y\left(0\right)$. When there is no noise, the equations are:

$$\frac{dp}{dt}=yp,\frac{dy}{dt}=-\gamma y.$$

(268)

Then the price is given by the expression,

$$p\left(t\right)=p\left(0\right)e^{y\left(0\right)/\gamma }exp\left\{-\frac{y\left(0\right)}{\gamma }{e}^{-\gamma t}\right\}.$$

(269)

With time, the drift term y tends to zero,

$$y\left(t\right)=y\left(0\right)e^{-\gamma t}\to 0(t\to \infty ),$$

(270)

and the price tends to the value,

$$p\left(t\right)\simeq p\left(0\right)e^{y\left(0\right)/\gamma }(t\to \infty ),$$

(271)

named the fundamental price. Thus the market price, it seems, should tend with time to equilibrium. Random noise induces only small fluctuations around the equilibrium fundamental price value.

Recall that this conclusion is based on the assumption of market efficiency. However, the questions remain: Are markets efficient in the sense discussed above? Are they in an equilibrium state or at least tend to equilibrium? Are they sufficiently stable? Are typical agents, acting in markets, really rational?

4.3. Herding Market Model

It has been noticed long time ago that agents are not completely rational, but their rationality is bounded [174]. The agents are always subject to hesitations, biases, superstitions, and other feelings [175]. The agents cannot be absolutely rational, since:

• they have numerous cognitive biases;
• they do not possess all necessary information;
• they are not able to accomplish instantaneous calculations;
• atually, they do not maximise some expected utility or other functionals.

Generally speaking, all agents, one way or another, are subject to various emotions that can essentially influence their decisions. It is hardly possible to numerically evaluate the role of emotions for each particular agent, since decision making is actually a random process [176], but the quantification of emotions during this process is admissible on aggregate level [177].

In particular, the members of a society are strongly subject to the so-called herding effect. As Poincaré [178] has written, “Individuals who are close to each other, as they are in a market, do not take independent decisions —they watch each other and herd.” The herd instinct acts similarly to feedback field in nonlinear systems, which in a financial market can lead to strong fluctuations [179,180,181] and even can trigger business cycles [182,183,184]. In addition, a market is not compulsorily equilibrium due to government interference, different external manipulations, corruption, insider trading etc. Financial crises, accompanied by strong fluctuations, as booms and crashes, can be caused by herding effect [28].

The behaviour of markets cannot be always stationary. Actually, sometimes markets look as almost stationary and efficient, but sometimes they become strongly nonequilibrium and inefficient. The market boom-crash fluctuations remind us the so-called heterophase fluctuations in statistical systems, where these fluctuations are caused by collective effects [185,186]. To make the market behaviour richer, so that to include in its description regime switching between conventions and business cycles, it is necessary to take account of collective effects [187].

Let us define the logarithmic mispricing,

$$x\equiv log\frac{p}{p^{*}},$$

(272)

where the price, p, is normalised to a fundamental price, $p^{*}$, and the logarithm can be taken with respect to any convenient base. The equation for the mispricing,

$$dx=ydt+{\sigma }_{1}dW$$

(273)

follows from the diffusion price model, as in the previous section. However the drift term is assumed to satisfy the stochastic equation,

$$dy=fdt+{\sigma }_{2}dW,$$

(274)

with its own drift force $f=f(x,y)$. As before, W is a Wiener process.

The drift force, $f(x,y)$, can be modelled by expanding $f(x,y)$ in powers of the variables x and y and then resorting to exponential summation [143], taking into account that the force, $f(x,y)$, has to be antisymmetric with respect to the replacement $x\mapsto -x$ and $y\mapsto -y$. The resulting drift force takes the form,

$$f(x,y)=\alpha x+\beta y+A{x}^{3}exp\left(-\frac{x^2}{{\mu }^2}\right)+B{y}^{3}exp\left(-\frac{y^2}{{\lambda }^2}\right).$$

(275)

The meaning of the terms in the force (275) is as follows.

The term $\alpha x$ describes the individual response to changing price, with $\alpha <0$ being correcting response, while $\alpha >0$ being speculative response. The term $\beta y$ characterises the individual response to the changing price drift, with $\beta <0$ being the correcting response named market friction.

The term

$$A{x}^{3}exp\left(-\frac{x^2}{{\mu }^2}\right)$$

describes collective response to changing price, with $A<0$ corresponding to the correcting response and $A>0$ to the speculative response. The parameter $\mu$ defines the measure of uncertainty in the price value. The absence of uncertainty implies a fully informed market, when $\mu \to 0$, hence

$$A{x}^{3}exp\left(-\frac{x^2}{{\mu }^2}\right)⟶0(\mu \to 0).$$

(276)

The herding behaviour corresponds to $A>0$.

The collective response to the varying price drift is modelled by the term

$$B{y}^{3}exp\left(-\frac{y^2}{{\lambda }^2}\right).$$

Here, $B<0$ corresponds to the contrarian response, while $B>0$ corresponds to the speculative response. The parameter $\lambda$ measures the level of market freedom, or deregulation. If a market is over-regulated, with $\lambda \to 0$, then there is no collective response,

$$B{y}^{3}exp\left(-\frac{y^2}{{\lambda }^2}\right)⟶0(\lambda \to 0).$$

(277)

Herding occurs under $B>0$.

The phase portrait, defining $y=y\left(x\right)$, is given by the equation

$$\frac{dy}{dt}=\frac{f(x,y)}{y}.$$

(278)

Depending on the relative value of the model parameters, there exists a rich variety of different regimes describing markets with coexisting equilibrium, conventions, and business cycles. The inclusion of noise defines stochastic dynamics that is characterised by nonlinear geometric random walk processes with spontaneous regime shifts between different conventions and business cycles. This model provides a natural framework to explain dynamic regime shifts between different market states. These shifts are the result of the interplay between the individual and herding effects, under the presence of noise [187]. The change of dynamic regimes induced by external noise leads to noise-induced dynamic transitions [188].

4.4. Time Series Analysis

Financial and economic time series are a particular kind of time series, whose analysis is necessary in many applications; the standard methods of the series analysis can be found in Refs. [189,190,191]. In this Section, a method that has been developed rather recently [192] is presented. This method is based on self-similar approximation theory (see reviews [193,194]) interpreting the transfer between approximants in approximation space as the motion of a dynamical system, with the approximant order playing the role of discrete time. The idea of the approach is the assumption that, in given data, filtered from noise, there exists a hidden law of evolution, which gives the possibility of understanding the future system development with a finite forecast horizon.

Suppose one observes a series of financial or economic market data expressed in the form of numbers, e.g., values or market indices, resulting in a set of data $z_i$ at the related moments of time $t_i$. It is necessary to keep in mind that the raw data as such are not representative for characterising the market dynamics, since these data may contain a great deal of noise. Therefore there is no clear convincing reason to analyse time series data in precisely the form in which they are provided [191]. The existence of noise is typical for practically all large systems, whether in physics or economics. Even equilibrium systems can exhibit rather strong fluctuations (see, e.g., [195]).

In markets, generally, there can happen two kinds of noise, caused by exogenous or endogenous sources. The former can be due to asset and wealth shocks caused by wars or other disasters, abrupt population changes, and such [196]. Endogenous noise is produced by the system itself. The evolution of an economic system is essentially a disturbance of equilibrium in the economy. Within each economic system there always exists a source of energy sufficient for disrupting any equilibrium [197]. A long-term equilibrium cannot be reached, but fluctuations and noise are everlasting [198,199]. There exist different ways of filtering out noise [189,190,191]. The simplest and rather reasonable way is as follows.

Step 1. Let us separate the studied time period into several time intervals, e.g., $k+1$ intervals, enumerated by the index $n=0,1,2,\dots ,k$, with the corresponding market data sets,

$$\{z_i^{\left(n\right)},t_i^{\left(n\right)}:i=1,2,\dots M_n\}(n=0,1,2,\dots ,k).$$

(279)

Depending on the scenario under study, each interval can be a day or a week or what is more appropriate.

The values $z_i^{\left(n\right)}$ are subject to random noise. This noise can be smoothed out by composing, for each time interval n, the average values,

$$x_n\equiv \frac{1}{M_n}\sum _{i=1}^{M_n}{z}_i^{\left(n\right)}.$$

(280)

The value $x_n$ can be ascribed to a point of time $t_n$ inside the nth time interval, which, e.g., can be chosen as the average of $t_i^{\left(n\right)}$, or as the left or right time point of the time interval in Equation (279), or, if the time intervals are equal to $\Delta t$, as $t_n=n\Delta t$. Then, one obtains the smoothed data set of the averages,

$${\mathbb{D}}_k=\{x_n,t_n:n=0,1,2,\dots k\}.$$

(281)

Step 2. Assume that the data set (281) represents a function $f\left(t\right)$ passing through the given time points. The most complete representation of $f\left(t\right)$ is given by the polynomial [200,201],

$$f_k\left(t\right)=\sum _{n=0}^k{a}_n{t}^n=a_0\left(1+\sum _{n=1}^k\frac{a_n}{a_0}{t}^n\right),$$

(282)

such that

$$f_k\left(t_n\right)=x_n(n=0,1,2,\dots ,k),$$

(283)

which defines $a_n=a_n(x_1,x_2,\dots ,x_k)$.

Step 3. The polynomial (282) can be extrapolated beyond the time $t_k$ by employing self-similar approximation theory [193,194], using self-similar exponential approximants [142,143]. Then, one obtains the approximants of first order,

$$f_1^{*}\left(t\right)=a_{0}exp\left(c_{1}t\right),$$

(284)

where $c_1$ is defined by comparing the expansion in powers of t of Equation (284) with the first order series (282), which gives

$$c_1=\frac{a_1}{a_0}.$$

(285)

The second-order approximant reads:

$$f_2^{*}\left(t\right)=a_{0}exp(c_{1}texp\left(c_{2}t{\tau }_2\right)),$$

(286)

with the control function ${\tau }_2$ defined by the fixed-point equation,

$${\tau }_2=exp\left(c_{2}t{\tau }_2\right),$$

(287)

and with $c_2$ following from the accuracy-through-order procedure by equating the expansion of the approximant (286), under ${\tau }_2=1$, with the polynomial (282), which leads to

$$c_2=\frac{a_2}{a_1}-\frac{a_1}{2a_0}.$$

(288)

The third-order approximant is

$$f_3^{*}\left(t\right)=a_{0}exp(c_{1}texp(c_{2}texp\left(c_{3}t{\tau }_3\right))),$$

(289)

with the parameter $c_3$ defined by the accuracy-trough order procedure,

$$c_3=\frac{a_3}{a_1{c}_2}-\frac{c_2}{2}-c-\frac{c_1^2}{6c_2},$$

(290)

and the control function ${\tau }_3$ prescribed by the fixed-point condition,

$${\tau }_3=exp\left(c_{3}t{\tau }_3\right).$$

(291)

Continuing this procedure, one gets:

$$f_k^{*}\left(t\right)=a_{0}exp(c_{1}texp(c_{2}t\dots exp\left(c_{k}t{\tau }_k\right))\dots ),$$

(292)

with the control functions

$${\tau }_k=exp\left(c_{k}t{\tau }_k\right).$$

(293)

Generally, control functions, ${\tau }_k={\tau }_k\left(t\right)$, can be defined by fixed-point conditions [202], as above, or by minimising a cost functional [203].

Extrapolating the approximants (292) to time $t>t_k$ gives nonlinear forecasts that can describe many booms and crashes which cannot be explained by standard methods; see detailed description of a number of examples in Refs. [202,204,205,206].

It is necessary to keep in mind that market bubbles and crashes can be provoked by many causes. There are indefinitely large number of functional relationships that can lead to such unstable effects as bubbles and crashes [28,207]. The standard time-series autoregressive methods are based on presumed structural stability and linearity, because of which they cannot in general predict unstable nonlinear, and strongly nonequilibrium market movements, such as booms and crashes [208]. Markets are largely influenced by human emotions resulting in herding behaviour [208,209,210]. In the words of Woods [210], “It is easier to forecast weather than to predict stock market prices”, and the essence of the standard autoregressive models is “garbage in, garbage out”.

The exponential extrapolation, described above, is highly nonlinear and gives hope for the possibility of grasping nonequilibrium market dynamics. However, since the origin and nature of different booms and crashes can be rather different, it may happen that it is necessary to employ different variants of self-similar extrapolation, for instance, different conditions defining control functions.

5. Conclusions

The content of this paper was based on the lectures that have been given by the author during several years at the Swiss Federal Institute of Technology in Zürich (ETH Zürich). The first part [4] considered equilibrium social systems. In this part of the lectures, nonequilibrium social systems are discussed. In addition to the general information, several novel evolution equations were studied, describing phenomena such as punctuated evolution, symbiosis of species with interactions through functional carrying capacity, dynamic phase transitions caused by the variation of growth rates, the stability of self-organised societies, and herding in markets.

The whole theme of social physics is huge and, certainly, it is impossible to cover it in a single survey. The reader can find many more models in the cited literature. Here, the choice of the touched-on topics was motivated by the interests of the author.