Selected Topics of Social Physics: Nonequilibrium Systems

This review article is the second part of the project ``Selected Topics of Social Physics". The first part has been devoted to equilibrium systems. The present part considers nonequilibrium systems. The style of the paper combines the features of a tutorial and a review, which, from one side, makes it easy 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. The present material is based on the lectures that the author had been giving during several years at the Swiss Federal Institute of Technology in Zurich (ETH Zurich).


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 researches.
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 characterized by evolution equations describing how the systems develop and what are their evolutionary stable states.
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 behavior of social systems. Second, to illustrate these methods on 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 characterizing symbiosis is presented. The peculiarity of constructing the evolution equations for structural self-organized societies is discussed. Some examples of financial markets are treated, emphasizing the important role of herding behavior.

Dynamical Social Systems
The evolution of nonequilibrium social systems is described by the time dependence of observable quantities. The time dependence is usually given through differential equations. Sometimes, one considers discrete time, when the evolution equations are given by difference equa-tions. Generally, time is a continuous variable. Here we shall mainly deal with differential equations.

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 characterizing the system evolution, from the mathematical point of view, is a dynamical system which many books are devoted to (see, e.g., [104][105][106][107][108]). The main definitions we shall need are briefly surveyed below.
Let the set of observable quantities of interest be x 1 , x 2 , . . . , x N . Each of the observables varies inside the related manifold X n , which is a part of the real-number axis R = (−∞, +∞). In order not to overload the text by notations, we shall simply write that x ∈ R N .
The set x ≡ {x n : n = 1, 2, . . . , N} ∈ R N (2.1) is the system state. The manifold 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 usually admissible to set t 0 = 0. So, often time is treated as the variable in the interval [0, ∞). 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 ) , with a state x 0 given at a moment of time t 0 , called initial state, 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 ′ , x 0 ) = x(t, x(t ′ , x 0 )) . The states at fixed times are called the points of the trajectory. General dynamical systems in continuous time are usually represented by differential equations in the normal form where f (x, t) = {f n (x, t) : n = 1, 2, . . . , N} . (2.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 (2.7) is termed a nonautonomous equation, while when f (x, t) = f (x) does not depend on time, this is an autonomous equation. The intermediate case is when at finite times the evolution equation is nonautonomous, but, as time increases, the equation becomes autonomous, so that Then the equation is called transiently nonautonomous. 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, . . . , N} , (2.9) for which there exists a time t * , such that The norm of a dynamic state is denoted as which is termed vector norm or Euclidean norm. This norm describes the character of motion that can be singular, or collapsing, and it can be nonsingular, or noncollapsing. The motion, starting at a point x 0 , is strongly singular, if there exists such a finite moment of time t c = t c (x 0 ), called critical, when the trajectory diverges to infinity, The motion is weakly singular, when the trajectory diverges as time goes to infinity: The motion is nonsingular, if the trajectory is bounded and never diverges.

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, Global Lagrange Stability. The motion is globally Lagrange stable, if it is Lagrange stable for all initial conditions, so that The solution x(t, x 0 ) is Poincaré stable if, for any ε > 0, there exist δ ε and a time t ε , called recurrence time, such that, Since any moment of time can be accepted as a new initial point, there can exist an infinite sequence or recurrence times.
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 ∈ R N . Asymptotic Lyapunov Stability. The solution x(t, x 0 ) is asymptotically stable, if it is Lyapunov stable and there exists δ ε > 0 such that if 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 ∈ R N .
Definitions of stability can be straightforwardly reformulated for fixed points. For the autonomous equations, fixed points x * are given by the condition Lyapunov Stable Fixed Point. The fixed point x * is Lyapunov stable, if for any ε > 0 there exists δ ε > 0 such that, if x 0 − x * < δ ε , then for all t > 0.
Globally Asymptotically Stable Fixed Point. The fixed point x * is globally asymptotically stable if it is asymptotically stable for all initial conditions x 0 ∈ 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.

Method of Linearization
Lyapunov developed two methods for analyzing the stability of motion. The most often used is the linearization method that we briefly remind below.
The analysis of stability starts with considering a small deviation δx(t) from the evolution of the variable x(t): with the Jacobian matrixĴ We can define the eigenvalues and eigenvectors of the Jacobian matrix by the eigenproblem The linearization is usually done 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, Then, the Jacobian, as well as its eigenvalues and eigenvectors, in the limit of large time, do not depend on time:Ĵ from a fixed point, given by the definition The Jacobian eigenproblem takes the form The solution to the linear equation for the deviation can be represented as an expansion over the eigenvectors of the Jacobian matrix: The Lyapunov exponents, or characteristic exponents, related to a fixed point, are the real parts of the Jacobian eigenvalues evaluated at this fixed point, Lyapunov exponents characterize the asymptotic stability of fixed points and describe the type of the asymptotic dynamics. The sign of the largest Lyapunov exponent (2.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. But 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 literature [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 For each fixed point, there is a single Lyapunov exponent The fixed point is asymptotically stable, if λ < 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 λ = 0. In this situation, the fixed point can be a center, around which the trajectory oscillates. The fixed point is asymptotically unstable, when λ > 0. Then, there is no basin of attraction for this fixed point, but all trajectories diverge from it.

Plane Motion
The motion on a plane is described by a two-dimensional dynamical system. That is, the phase space is R 2 . For an autonomous system, there are two differential equations with initial conditions Dividing the second equation over the first gives the relation defining the trajectory y = y(x, x 0 , y 0 ) .
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 ) : The Jacobian is a two-by-two matrix with the elements

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 characterized by differential equations. Let us consider a society composed of several different groups of populations enumerated by the index i = 1, 2, . . .. In an i-th group there are N i members. The number of the members N i = N i (t) 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 In other cases, N can represent a company capitalization, with N i being the capitalization of the company parts. The number of the society members can be very large because of which it is more convenient to use reduced numbers, for example normalizing N i by a fixed scaling N 0 , so that Often N 0 can be chosen as the total number of the society members at the initial time, N(0), 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, . . . of all variables forms the society state. The evolution equations for the society groups usually have the form of the rate equation is the effective rate of the population group change and Φ i = Φ i (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 where γ 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 as and in the reduced form the evolution equation becomes with the notation

Delay Differential Equations
Sometimes, modeling 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 with different discrete delays τ n for each population, reads as where τ n are delay times and f = {f n } is also a vector. The history conditions for negative times are to be given by prescribed history functions with h n (t) being explicit functions of time.
If for each species, the delay time τ n = τ is the same, the evolution equation takes the form The solution of a delay equation can be done following the step-by-step method. For instance, let us consider the case of a single population described by equation (2.55), with the history condition x(t) = h(t) for t ≤ 0. One considers the solution separately in the subsequent time intervals [0, τ ], then [τ, 2τ ], then [2τ, 3τ ], etc. At the first step, the solution is defined by the equation with the initial condition At the second step, denoting one finds the solution from the equation with the initial condition Continuing this procedure for subsequent time intervals, at the step k, the solution is given by the equation with the initial condition In practice, the solution at each step is calculated numerically.

Examples of Evolution Equations
There exist many variants of evolution equations. Below, we briefly mention several the most known equations.

Malthus Equation
Malthus [116] was interested in the Earth population growth. He noticed that the Earth population varies according to the law The variation rate is the difference between the growth, γ + , and death, γ − , rates of the population, γ = γ + − γ − . As is evident, with time, the population either diminishes to zero, if γ < 0, or tends to infinity, if γ > 0. Since for the people on the Earth, the birth rate surpasses the death rate, hence γ > 0, the population will grow to infinity. The singular behavior 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 γ equal to zero.

Neo-Malthusian Catastrophe
Moreover, some authors have suggested that the situation is even worse than predicted by Malthus, and the Earth population can become strongly singular, diverging at a finite moment of time. This situation 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 Then the evolution equation reads as With the scale N 0 = N(0), we have the reduced equation The solution shows that the population diverges as

Logistic Equation
Verhulst [118] 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 taking into account the existence of the carrying capacity K entering the term characterizing the competition between the society members. The Malthus equation can be valid only at the very beginning of the population growth when N ≪ K.
The logistic equation has found a number of applications describing the growth of population in demography, company capitalization in economics, neuron signals in neural networks, tumor size in medicine, reactant mass in chemistry, and so on. The reduced form of the equation, with the population normalized to N 0 = N(0), is The solution to the logistic equation is the logistic function also called sigmoid function. Its limit with respect to time is finite, hence there is no singular behavior. From the point of view of the Lyapunov analysis, fixed points are defined by the equation resulting in two fixed points, x * 1 = 0 and x * 2 = γ/a. The Jacobian The model [119,120] 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 where all parameters are assumed to be positive. As usual, the initial conditions are denoted as When x 0 = 0, then y(t) dies down to zero as time increases. If y 0 = 0, then x(t) explodes to infinity with time. More interesting are the nonzero initial conditions. Let us employ the Lyapunov stability analysis. The fixed-point equations yield two fixed points: one is trivial 83) and the other is nontrivial, The Jacobian matrix is composed of the elements which lead to the Jacobian trace and determinant At the trivial fixed point, we have The Jacobian eigenvalues are Hence, the trivial fixed point is unstable, being a saddle point. For the second fixed point, we get hence the Jacobian eigenvalues become This means that the nontrivial fixed point is an elliptic point for any nonzero parameters of the evolution equation. The corresponding temporal behavior is given by permanent oscillations. The predators and prey oscillate with the same frequency ω, 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.

Singular Malthus Equation
The Malthus equation (2.65) can be modified [121] to the form whose solution diverges by power law if ε = m − 1 > 0 and the critical time is For m = 1, we return to the original Malthus equation that yields an exponentially divergent solution at t → ∞. Such strongly singular solutions were applied to rationalize the super-exponential growth of the human world population, financial markets, material failures, earthquakes, climate changes, and dynamics of other systems (see [122]).

Generalized Lotka-Volterra Model
The Lotka-Volterra model can also be generalized to the case of multiple species for the populations N i , with i = 1, 2, . . ., yielding The dimensionless populations can be introduced by normalizing N i , for example, to the total initial population N(0) = i N i (0), which gives the reduced form The signs of the parameters γ i and a ij can be different, providing a large variety of possible solutions [123]. For two species, we return to the original Lotka-Volterra model.

Predator-Prey Kolmogorov Model
A particular form generalizing the Lotka-Volterra Model, in which is called the predator-prey Kolmogorov model [124,125]. This model can again be rewritten in the reduced form To model concrete cases, it is required to specify the functions f 1 and f 2 .

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 , 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 (N 2 ). The supply function is taken in the form, such that the nutrient becomes depleted (N 2 → 0), as time increases, while the bacteria population reaches a fixed point characterizing a stationary value [126].

Hutchinson Delayed Equation
There are as well several generalizations of the logistic equation (2.76). Thus [127] considered an effective reproductive rate delayed in time, which gives the equation Here K is a fixed carrying capacity and τ is a delay time. The solution to this equation displays oscillations above the logistic curve.

Multiple-Processes Delayed Equation
There exist many variants of the delay logistic equations [113,128] designed for the description of single-species population N. For example, one can introduce multiple carrying capacities K j and multiple delay times τ j characterizing different processes in the dynamics of a population. The single-species population dynamics, with multiple processes reads as

Peschel-Mende Hyperlogistic Equation
The accelerated growth of population, such that exists for the human world population, can be described [129] by means of two additional powers m and n, which leads to the equation 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 exceed the given carrying capacity K. The Verhulst logistic equation (2.76) is recovered for m = 1 and n = 1.

Hyperlogistic Delayed Equations
The Peschel-Mende hyperlogistic equation can be extended to the delayed equation [130] This equation can also be treated as an extension of the Hutchinson delayed equation (2.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.

Replicator Equation
A very popular equation employed for modeling 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 characterizes 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 is the average society fitness. The total population is assumed to be constant, hence the equation is defined on a simplex [123].
In terms of the reduced fractions 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 the 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, say There are two fixed points, x * 1 = 0 and x * 1 = 1. The Jacobian is At the fixed points, this gives 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 However, when the fitnesses are such that f 1 is smaller than f 2 , then only the second species survives, since This result one often formulates as a sentence "the strongest survives".

Free Replicator Equation
It is interesting that the same qualitative result follows from the replicator equation that is free from the normalization conditions (2.104) or (2.105), which can be called free replicator equation. Then we need to deal with the number of equations equal to the number of species. Thus, for the case of two species, we have two equations We assume that the species fitnesses are different, f 1 = f 2 , otherwise the species could not be distinguished.
There are three fixed points, where either At the fixed point (2.116), for the Jacobian we get thence the Laypunov exponents are This tells us that the fixed point (2.116) is stable provided that f 1 > f 2 , when the second species dies out, while the first species survives.

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 realistic situations, it may happen that the evolution of a society is subject to random perturbations, termed random noise. Then, instead of Eq. (2.7), one comes to the equation containing, in addition to f (x, t), a term describing the action of an external noise. Usually, setting since the noise term is often modeled by the Gaussian white noise that is not differentiable. In the case of white noise, the variable W (t) characterizes a standard Wiener process. For time in the interval [0, t], the random variable W (t) is drawn from the normal law with zero mean and the standard deviation √ t. The equation (2.126) is called stochastic differential equation (see the books [131][132][133][134]). The quantity f (x, t) is called drift term and g(x, t), diffusion term. For a while, we consider a single variable x = x(t). The generalization to many variables will be given at the end of the section.
In practical calculations, the dynamics of x(t), in the presence of noise, is usually represented through the Euler-Maruyama scheme [135]. For this purpose, the interval [0, t] is partitioned in discrete points t n = n∆t, with n = 0, 1, 2, . . .. Then equation (2.126) is integrated between t n and t n+1 = t n + ∆t, which gives (2.127) The first integral in the right-hand side is a usual Riemann integral that, for small ∆t, can be written as The Riemann integral does not depend on the choice of the variable t inside the interval [t n , t n + ∆t], so that instead of t n in the right-hand side of the above formula it is possible to take any t ′ n ∈ [t n , t n + ∆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 (2.127), following Stratonovich or Ito [131][132][133][134]). According to Stratonovich, for infinitesimally small ∆t, one takes [g(x(t n + ∆), t n + ∆) + g(x(t n ), t n )]/2, while according to Ito one has to take g(x(t n ), t n ). When employing the Euler-Maruyama numerical scheme [135] one uses the Ito representation. Then the second integral in (2.127), understood in the sense of Ito, because of which it is called the Ito integral, writes as (2.129) The random variable W (t n ) is drawn from the normal distribution with zero mean and standard deviation √ ∆t. For each t 1 < t 2 , the normal random variable W (t 2 ) − W (t 1 ) is independent of the random variable W (t 1 ). The difference can be represented as where R(t n ) is drawn from the normal distribution with zero mean and standard deviation one.
Thus we come to the iterative equation In those cases, where the diffusion term is constant, g(x, t) = σ, we have Under the situation with several variables x i , where i = 1, 2, . . ., and several noise terms, one has For a diagonal diffusion matrix g ij (x, t) = σ i δ ij , this reduces to The existence of noise superimposes on a dynamical trajectory random deviations, whose amplitude depends on the value of σ.

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 realization of the random noise. For the stochastic differential equation (2.126), the probability density of x is described by the Fokker-Planck equation [136] with the diffusion coefficient Introducing the probability current yields the continuity equation 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 (x), g(x), and D(x). Respectively, the probability density tends to p(x), satisfying the equation in which the stationary probability current is The latter equation implies that Suppose that the quantity x varies in an interval starting at x 0 , so that x ≥ x 0 . Assuming the reflecting boundary condition which is equivalent to the equation The solution to the latter equation is with the constant C defined from the normalization condition where the integration is over the whole range of variation of x. Expression (2.145) is named potential solution, as far as it can be rewritten in the form with the potential All above formulas can be easily generalized to a multivariate case, where the Fokker-Planck equation reads as with the diffusion matrix of the elements (2.150)

Generalized Evolution Equations
Usually, in the evolution equations the rates of change and population interactions are treated as constant parameters. In more realistic, and hence more general, situations they may be functions of the population variables. Below, we consider several such generalized equations.

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 where γ i is an effective birth rate for a population, if γ i is positive, or a death rate, when γ i is negative. If one considers not population evolution, but dynamics of production of a firm or like that, then a positive γ i takes the sense of gain rate, while a negative γ i means loss rate. The quantity A ij describes effective interactions between the i-th and j-th populations. Usually, A ij is treated as a constant or, sometimes, as a given function fixed by external forces [137][138][139]. Generally, in the population interactions it is admissible to distinguish direct interactions and indirect interactions through the mutual influence on their carrying capacities [122,140]. Therefore, the effective interactions can be represented in the form with B ij being a direct interaction parameter and K j , carrying capacity of the j-th population. The influence of populations on the carrying capacity is assumed, which implies that it depends on these populations: 3) The populations entering K j are delayed, since to induce a change in the carrying capacity requires time. The term Φ i is caused by the influx of populations from outside, if any. As always, it is more convenient to deal with reduced dimensionless quantities, for which we introduce Here the normalization quantity N 0 , for a while, is not defined. It will be chosen later on so that to simplify the reduced equation. Then the evolution equation reads as If there are no external fluxes from outside, the term ϕ 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 (3.5) reduces to Using the sign function Defining the dimensionless carrying capacity and measuring time in units of γ −1 yields the evolution equation Formally, this looks as a logistic equation. However the principal difference here is the functional dependence of the carrying capacity on the delayed population, q = q(x(t − τ )) . (3.10) 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

This is equivalent to the expansion
The carrying capacity has to retain its sense, hence to be finite for asymptotically small amount of population. This means that To simplify the equation, we may choose the normalization quantity N 0 in the form hence q 0 = 1, as a result of which the expansion of the carrying capacity becomes 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 [140]. 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 (3.11). For this purpose, we can resort to the exponential summation guaranteeing a positive effective limit [141][142][143]. This gives q(x) = exp{bx(t − τ )} . (3.15) In this way, the evolution equation reduces to [144] dx dt The production parameter b characterizes 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 situation, it is, actually, a destruction parameter. The initial condition for the delay equation is Depending on the signs of the parameters γ and B, there can happen the following four different types of evolution models characterized by:

Evolutionary Stable States
Delay equations allow for the Lyapunov stability analysis, similar to that for the standard differential equations [128]. The fixed points, or stationary states, are defined by the equation There always exists the trivial solution existing for any sgnγ and sgnB. Nontrivial solutions require the validity of the relation which imposes the constraint Under this condition, fixed points are given by the equation When the population destroys its carrying capacity, that is b < 0, there can exist one fixed point in the range 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 3.23) and the other for (3.24) At the bifurcation point b = 1/e, the fixed points x * 2 and x * 3 coincide: There are no stationary nontrivial states for b > 1/e. The fixed-point stability is characterized by the behavior of small deviations from the fixed point. Substituting into the evolution equation (3.16) the definition yields the linearized equation When the solution to (3.27) is bounded, the solution to equation (3.16) is Lyapunov stable. When the solution to (3.27) converges to zero for t → ∞, a fixed point is asymptotically stable.
The trivial fixed point x * 1 = 0 is stable when but it is always unstable for sgnγ = 1, The motion in the vicinity of nontrivial fixed points is described by equation (3.27), which, using relation (3.19), reduces to the equation Looking for the solution in the form we get the equation for the Lyapunov exponent λ depending on the studied fixed point.

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 [145][146][147]. 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 [144]. This model might characterize not only the evolution of biological species, but also the evolution of firms and other organizations.
In order to emphasize 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 characterized by gain and competition [144], where sgnγ = sgnB = 1 .
In that case, the trivial fixed point x * 1 = 0 is never stable. The nontrivial fixed point x * 2 is stable when either where 34) or in the range 1 e < x * 2 < e −e < b ≤ 1 e , τ > 0 . The fixed point x * 3 > e is always unstable under gain and competition. So that the sole stable fixed point is x * 2 ∈ (0, e). When the production parameter b is seminegative, 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 [144].
Theorem on population boundedness. The solution x(t) to the evolution equation (3.32), for b ≤ 0, any finite τ ≥ 0, and any initial conditions x 0 ≥ 0, is bounded for all times t ≥ 0, and, for b < 0, there exists a time t 0 = t 0 (x 0 , τ ) such that The proof is given in Ref. [144].
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 semiline x 0 ≥ 0, but a limited interval. This happens for b ∈ (0, 1/e), when the basin of attraction is given by the inequalities For the values of b satisfying condition (3.36), but with x 0 > x * 3 , the solution x tends to infinity, as t → ∞.
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 > 1 e , τ > 0 , x 0 > 0 . 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 The unbounded punctuated behavior 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 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 there appear sharp reversals after almost horizontal plateaus. For b < −e , τ ≤ τ 2 , x 0 > 0 , (3.41) 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 , τ ≥ τ 2 , x 0 > 0 , (3.42) then there develops the regime of everlasting oscillations.
The discussed regimes of dynamics under gain and competition are described in detail in Ref. [144], 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 organizations, growth of alive organisms, improvement of individual abilities, proliferation of cells and bacteria, as well as in the decay of biological organisms and societies [122,140,144,145].

Symbiosis of Species
The term symbiosis characterizes 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 in the books [148][149][150][151][152]. 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.

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 normalize the populations N i with respect to different normalization constants M i , introducing the fractions Then the reduced form of the equations, without external forces, can be written as The convenience of using different normalizations 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 [153]. Essentially different numbers of different populations can require different normalization constants.
The coexisting species are characterized 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, (3.45) Therefore equation (3.44) reduces to It is possible to introduce dimensionless carrying capacities that are functions q i = q i (x 1 , x 2 , . . .) such that q(0, 0, . . .) = 1, which is possible by defining the appropriate normalization constants K 0 i . Generally, the variables entering q i should be delayed in time, such that x i = x i (t − τ i ).
Let us consider the case when the parameters γ i and B ii are positive. And let us choose the normalization factors M i as Considering the case of two species, we define x 1 ≡ x and x 2 ≡ y. Then equations (3.49) reduce to the system of two equations Let us measure time in units of 1/γ 2 and define the ratio Thus we come to the equations (3.52) Now we need to define the carrying capacities q i that, generally, are the functions of the delayed populations For the case of two populations, q i = q i (x, y). If we expand the function q i (x, y) in a series in powers of populations, we get the form as 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 [122]. 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 strong mutual influence, the expansions can be summed by means of exponential summation [141][142][143], thus avoiding spurious zeroes in the effective carrying capacity [154,155]. In that way, we can derive the form (3.55) 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 characterizes the impact of the populations on the carrying capacity, when the symbiotic populations correlate with each other. From the general expression (3.55), 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 (uncorrelated symbiosis) , (3.56) 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 situation describes correlated symbiosis, when the species influence the carrying capacities being mutually correlated, which for two species reads as q 1 = e bxy , q 2 = e gxy (correlated symbiosis) . (3.57) 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 we give the classification of simbiotic types. In order to be precise and not to disturb the meaning, the formulation of the definitions below are given closely following Refs. [154,155].
(i) Mutualism is a relationship between different species when both of them receive mutual benefit, which corresponds to the symbiotic parameters 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 (3.59) Parasitism is a relation, when at least one of the coexisting species is harmful to the other one, which is characterized by one of the conditions As is seen, there exist various types of symbiosis described by the systems of differential equations that can be solved numerically [156]. Since symbiosis is extremely widespread in nature, there are numerous particular examples of coexisting species among biological and ecological societies, including bacteria and viruses [148][149][150][151][152]157].

Uncorrelated Symbiosis
Dynamics of symbiotic species, under the uncorrelated symbiosis and the assumption of approximately equal rates γ 1 and γ 2 , are described by the system of equations with the symbiotic parameters b ∈ (−∞, ∞) and g ∈ (−∞, ∞). For simplicity, we assume that the delay times of the populations entering the carrying capacities are very small, so that can be neglected. The initial conditions are x 0 = x(0) and y 0 = y(0).
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 that can be rewritten as The Lyapunov exponents are defined by the expressions 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 situation 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.
The details can be found in Ref. [155].

Correlated Symbiosis
Correlated symbiosis of two species is characterized by the system of equations dx dt = x − x 2 e −bxy , dy dt = y − y 2 e −gxy , (3.70) where we again assume that γ 1 is close to γ 2 . Similarly to the previous 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 bx * y * , y * = e gx * y * , (3.71) that can be rewritten as The Laypunov exponents are defined by the relations In the case of mutualism, where b > 0 and g > 0, there are the following types of behavior: (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 is positive, or both parameters are negative, the following situations can happen: (i) Convergence to stationary states.
(ii) Unbounded growth of parasitic population and dying out of host population.
The detailed investigation is given in Ref. [155].

Mixed Symbiosis
In the previous examples describing symbiosis, 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 nonsymmetric relations that can be characterized 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 Similarly to the previous cases, 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 that can be represented as The Lyapunov exponents are 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 situations: (i) Convergence to stationary states.
Details can be found in Ref. [155]. Applications to financial markets, considering symbiosis between asset prices and bonds, resulting in periodically growing and collapsing bubbles, are analyzed in Ref. [158].

Role of Growth Rates
When the growth rates γ 1 and γ 2 are essentially different, it is necessary to consider the system of equations in which α ≡ γ 1 /γ 2 and time is measured in units of γ −1 2 . Without the loss of generality, it is possible to take α > 1. In the case of uncorrelated symbiosis, we have There are three trivial fixed points: the unstable node {0, 0}, with the Lyapunov exponents λ 1 = 1 and λ 2 = α; a saddle {1, 0}, with the Lyapunov exponents λ 1 = 1 and λ 2 = −α; and the saddle {0, 1}, with the Lyapunov exponents λ 1 = −1 and λ 2 = α. The nontrivial fixed points do not depend on the value of α and are defined as above. The related Lyapunov exponents are Although the fixed points do not depend on the rate α, but the Lyapunov exponents do depend on it, as well as the related basins of attraction also depend on α, which hence influences the stability of fixed points [159]. Therefore by varying α it is admissible to obtain different regimes of motion, hence to realize dynamic phase transitions by the sole variation of the growth rate, with keeping other parameters unchanged.
A usual dynamic transition implies a qualitative change of dynamical behavior, 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 behavior 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 behavior 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 it [160].

Self-Organized Society
The evolution of biological societies, including human societies, is a principally important old problem studied in voluminous literature beginning with Darwin [161,162]. Societies usually are 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 usually studied on the basis of the replicator equation discussed in the above sections. As has been explained, 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 evident that for a closed self-organized 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.
Our aim is to consider a self-organized society, whose means of survival are produced inside the society itself. A closed self-organized 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.

Trait Groups
In the present section, a self-organized structured society is described composed of the trait groups representing four types of typical agents: cooperators, defectors, regulators, and outsiders. The principal novelty of the approach is that, instead of a single fitness or utility for each group, we introduce relative utilities for each group with respect to the society as a whole, and mutual utilities with respect to each other [163].
Let us consider a society composed of several groups, whose fractions are defined as the ratios with respect to a normalizing number N 0 that can be chosen to define the initial total number of the society members However, at the later times such a normalization, 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. [164,165].
• Cooperators (x 1 ), 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 (x 2 ), 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 (x 3 ), 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 (x 4 ), 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 (3.81), it is possible to accept the rate equations These equations are formally equivalent to the rate equations (2.52) or to the Lotka-Volterra equations (2.93). However, the meaning of the parameters entering equations (3.84) is different. In the Lotka-Volterra equations, γ i are the birth-death rates and a ij are intensities of interactions. In our case, each γ 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 j-th group with respect to the i-th 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. [164,165].
Cooperators. They 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. Summarizing this, we have: Defectors. They 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: Regulators. They 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: Outsiders. They are not useful for the society. But 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: Taking into account the signs of the parameters, equations (3.84) write as 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, we shall not complicate the situation and will treat a ij 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 we shall name a ij simply as utilities.
The cooperators are the sole group producing resources for the whole society. These resources are denoted through γ 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 | = γ 1 . (3.90) The diagonal elements |a ii | describe the strength of competition among the members of a group i for the available for them resources γ i playing the role of a carrying capacity. The standard dependence for the term characterizing competition on the carrying capacity is the inverse dependence where C is a constant. From relation (3.90) it follows that C = γ 2 1 . Therefore 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 [166,167] that can be written as Also, there should exist the reciprocity relation signifying that, if a group i receives some utility from a group j, then the group j looses the same amount of utility. We 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 95) and the dimensionless influx of outsiders We measure time in units of γ −1 1 . Then equations (3.89) acquire the form of the system of equations with the right-hand sides 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, (3.99) In dimensionless units this is equivalent to the inequality If this condition becomes invalid, this means that defectors, regulators, and outsiders consume more than cooperators produce. This could be realized only under a supply from outside or by the destruction of the given carrying capacity, which implies an unstable situation.

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-organized society. Let us see what is the situation in our case of a self-organized society.
Considering cooperators and defectors, we have the evolution equations Looking for fixed points satisfying conditions (3.100), we find the sole solution Thus in a stable self-organized society, the amount of defectors cannot be larger than around 10%. If this amount is essentially larger, the society is not stable.

Three Coexisting Groups
The standard situation is the coexistence of three groups, cooperators, defectors, and regulators, while outsiders are not numerous, hence can be neglected. Then the equations are There is the sole fixed point satisfying the conservation law (3.100), which is stable when This shows that the minimal fraction of cooperators is the maximal fraction of defectors is max a,b x * 2 ≈ 0.08 (a ≈ 0.48, b = 0) , (3.110) and the maximal fraction of regulators is Again we 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. [163], where the influence of noise is also studied.

Models of Financial Markets
Financial markets can be treated as complex social systems, because of which their behavior 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 examples.

Efficient Market Model
The old standing hypothesis, the functioning of financial markets is based on, is the efficient market hypothesis. It is generally accepted [168,169] that the idea of efficient market has been anticipated by Bachelier [170] 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 [171,172] 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. Summarizing, 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 maximization of expected utility.
Let p = p(t) 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 in which y is a drift rate of the price, σ, 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 dp(t) dt = R(t)p(t) , with a randomly fluctuating rate where ξ(t) = dW (t)/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 While in the presence of noise, the numerical solution is given by the iterative scheme

Diffusion Price Model
In the diffusion price model, the drift term y in Eq. (4.1) is treated not as a constant, but as a random variable, so that the market is described by the system of equations dp = (ydt + σ 1 dW 1 )p , where γ > 0 is called market friction. The initial conditions are p(0) and y(0). When there is no noise, the equations are dp dt = yp , dy dt = −γy . Then the price is given by the expression With time, the drift term y tends to zero, and the price tends to the value 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?

Herding Market Model
It has been noticed long time ago that agents are not completely rational, but their rationality is bounded [173]. The agents are always subject to hesitations, biases, superstitions, and other feelings [174]. 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.
• Actually, they do not maximize some expected utility or other functionals.
Generally speaking, all agents, to more or less extent, 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 in principle a random process [175], but the quantification of emotions during this process is admissible on aggregate level [176].
In particular, the members of a society are strongly subject to the so-called herding effect. As Poincaré [177] 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 [178][179][180] and even can trigger business cycles [181][182][183]. 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 behavior 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 [184,185]. To make the market behavior richer, so that to include in its description regime switching between conventions and business cycles, it is necessary to take account of collective effects [186].
Let us define the logarithmic mispricing where the price p is normalized to a fundamental price p * and the logarithm can be taken with respect to any convenient base. The equation for the mispricing dx = ydt + σ 1 dW (4.12) follows from the diffusion price model, as in the previous section. However the drift term is assumed to satisfy the stochastic equation with its own drift force f = f (x, y). As usual, W is a Wiener process.
The drift force f (x, y) can be modeled by expanding f (x, y) in powers of the variables x and y and then resorting to exponential summation [142], taking into account that the force f (x, y) has to be antisymmetric with respect to the replacement x → −x and y → −y. The resulting drift force takes the form (4.14) The meaning of the terms in force (4.14) is as follows.
The term αx describes the individual response to changing price, with α < 0 being correcting response, while α > 0 being speculative response. The term βy characterizes the individual response to the changing price drift, with β < 0 being correcting response named market friction.
The term describes collective response to changing price, with A < 0 corresponding to correcting response and A > 0, to speculative response. The parameter µ defines the measure of uncertainty in the price value. The absence of uncertainty implies a fully informed market, when µ → 0, hence The herding behavior corresponds to A > 0. Collective response to the varying price drift is modeled by the term Here B < 0 corresponds to contrarian response, while B > 0, to speculative response. The parameter λ measures the level of market freedom, or deregulation. If a market is over-regulated, with λ → 0, then there is no collective response, Herding occurs under B > 0. The phase portrait, defining y = y(x), is given by the equation (4.17) 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 characterized 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 [186]. The change of dynamic regimes induced by external noise leads to noise-induced dynamic transitions [187].

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 their analysis can be found in the books [188][189][190]. Here we present a method that has been developed rather recently [191]. This method is based on self-similar approximation theory (see reviews [192,193]) 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 we observe 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 characterizing the market dynamics, since these data may contain a great deal of noise. Therefore there is no clear convincing reason to analyze time series data in precisely the form in which they are provided [190]. 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. [194]).
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 like that [195]. 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 [196]. A long-term equilibrium cannot be reached, but fluctuations and noise are everlasting [197,198]. There exist different ways of filtering out noise [188][189][190]. The simplest and rather reasonable way is as follows.
Step 1. Let us separate the studied time period into several time intervals, say k + 1 intervals, enumerated by the index n = 0, 1, 2, . . . , k, with the corresponding market data sets, The value x n can be ascribed to a point of time t n inside the n-th time interval, which, e.g., can be chosen as the average of t Step 2. Assume that the data set (4.20) represents a function f (t) passing through the given time points. The most complete representation of f (t) is given by the polynomial [199,200] f k (t) = k n=0 a n t n = a 0 1 + k n=1 a n a 0 t n , (4.21) such that f k (t n ) = x n (n = 0, 1, 2, . . . , k) , (4.22) which defines a n = a n (x 1 , x 2 , . . . , x k ).
Step 3. The polynomial (4.21) can be extrapolated beyond the time t k by employing selfsimilar approximation theory [192,193], using self-similar exponential approximants [141,142]. Then we obtain the approximants of first order f * 1 (t) = a 0 exp(c 1 t) , (4.23) where c 1 is defined by comparing the expansion in powers of t of (4.23) with the first order series (4.21), which gives c 1 = a 1 a 0 . (4.24) The second-order approximant reads as f * 2 (t) = a 0 exp(c 1 t exp(c 2 tτ 2 )) , (4.25) with the control function τ 2 defined by the fixed-point equation and with c 2 following from the accuracy-through-order procedure by equating the expansion of approximant (4.25), under τ 2 = 1, with polynomial (4.21), which leads to The third-order approximant is f * 3 (t) = a 0 exp(c 1 t exp(c 2 t exp(c 3 tτ 3 ))) , with the parameter c 3 defined by the accuracy-trough order procedure, 29) and the control function τ 3 prescribed by the fixed-point condition  Generally, control functions τ k = τ k (t) can be defined by fixed-point conditions [201], as above, or by minimizing a cost functional [202].
Extrapolating the approximants (4.31) to time t > t k gives nonlinear forecasts that can describe many booms and crashes which cannot be explained by standard methods. Detailed description of a number of examples can be found in Refs. [201,[203][204][205].
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,206]. The standard time-series autoregressive methods are based on presumed structural stability and linearity, because of which they cannot in principle predict unstable nonlinear, and strongly nonequilibrium market movements, such as booms and crashes [207]. Markets are largely influenced by human emotions resulting in herding behavior [207][208][209]. In the words of Woods [209], "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 the 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.

Conclusion
The content of this article is 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 are studied describing such phenomena as punctuated evolution, symbiosis of species with interactions through functional carrying capacity, dynamic phase transitions caused by the variation of growth rates, stability of self-organized societies, and herding in markets.
The whole theme of social physics is huge and, of course, it is impossible to cover it in a single survey. The reader can find a lot more of models in the cited literature. Here the choice of the touched topics is motivated by the interests of the author.
Funding: This research received no external funding.