Modelling Population Dynamics of Social Protests in Time and Space: The Reaction-Diffusion Approach

Abstract

Understanding of the dynamics of riots, protests, and social unrest more generally is important in order to ensure a stable, sustainable development of various social groups, as well as the society as a whole. Mathematical models of social dynamics have been increasingly recognized as a powerful research tool to facilitate the progress in this field. However, the question as to what should be an adequate mathematical framework to describe the corresponding social processes is largely open. In particular, a great majority of the previous studies dealt with non-spatial or spatially implicit systems, but the literature dealing with spatial systems remains meagre. Meanwhile, in many cases, the dynamics of social protests has a clear spatial aspect. In this paper, we attempt to close this gap partially by considering a spatial extension of a few recently developed models of social protests. We show that even a straightforward spatial extension immediately bring new dynamical behaviours, in particular predicting a new scenario of the protests’ termination.

1. Introduction

Understanding of the dynamics of riots, protests, and social unrest more generally is important in order to ensure a stable, sustainable development of various social groups, as well as the society as a whole. To achieve this goal, mathematical modelling has been gaining growing recognition recently as a powerful research approach [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. Indeed, since replicated sociological experiments with large groups of people are rarely possible (and it is hardly possible at all to simulate social unrest under controlled conditions), capturing the complexity of the social dynamics through tractable experiments is not feasible. Mathematical modelling and computer simulations create a “virtual society” where hypotheses can be tested and different scenarios investigated in detail, thus providing a feasible alternative to the experiment.

A variety of modelling approaches and techniques has been developed and used such as the costs-benefits analysis [9], social network models [5,10], agent based models [3,7], and behavioural and epidemiological models [2,11,13,14], both in spatial and non-spatial systems. The models range from relatively simple, allowing only for some basic feature of the phenomenon [3], to more complicated ones that take into account more details such as, for instance, heterogeneity of social norms and behavioural responses [9,15]. The case studies include a few well known events such as the Russian revolution of 1905–1907 [1], the French riots of 2005 [13], the London riots of 2011 [9], the Arab Spring [16], and the Yellow Vest Movement in 2018–2019 [17], in particular aiming to describe the intensity and timing of the unrest and to analyse the processes and mechanisms involved.

Social protests are known to exhibit a variety of different scenarios. For instance, the total number of people involved in protests (to which we refer to as “crowd”) can change through time in different ways, e.g., staying approximately constant for some time and then decaying fast, or decreasing steadily but slowly. Furthermore, the social structure of the crowd can change significantly, e.g., due to the individual variability of people’s perception and expectations [18]. Whilst a comprehensive description of protests requires several state variables and a multivariate analysis (e.g., to account for the individual variability, cf. [15]), one quantity of immediate importance is the total number of people participating in the given event at the given time. Arguably, understanding how the total number may evolve with time provides crucial information needed to ensure security measures and to minimize possible damage and violence. An inspection of the data on the number of people participating in protests and riots available for several events of social unrest [9,13,19,20,21] immediately reveals that there are some generic properties of the dynamics. There usually exists the main peak, i.e., the maximum in the participants number (sometimes followed or preceded by secondary peak(s) of smaller magnitude) followed, in the course of time, by a tail where the number eventually decays to either zero or to a small background value. However, the rate of decay at the tail can be significantly different for different events. Whilst in many cases, the number of protesters decays fast, so that the unrest is effectively over in a matter of weeks [9,13], there have been several cases where the rate of decay is very slow. One such case is given by the Russian revolution of 1905–1907 [19,20], and another case is given by the Yellow Vest Movement [21]. As a result of the slow decay in the number of protesters, the duration of unrest is considerable (about two years for the Russian revolution, about one year for the Yellow Vest Movement, at the time when this paper was written). A question arises here as to whether the events with significantly different decay rates are different altogether (e.g., controlled by different factors) or they are the manifestation of the same social dynamical phenomenon, but occurring in a different parameter range. A related question is whether it might be possible to develop a “universal” mathematical model that would embrace the protests with a significantly different decay rate. Addressing this question is one of the goals of this paper.

Along with the changes in time, in many cases, riots and protests also have a clear spatial aspect. For instance, the timing and the intensity of protests can be significantly different at different locations. Experience gained from studies on dynamical systems of a different origin, e.g., in ecology [22], indicates that coupling between local protests (e.g., in the context of social dynamics, by the information exchange and/or spread of ideas) can organize them into a global spatiotemporal pattern. The simplest example of a spatiotemporal pattern is the travelling wave (sometimes also referred to as the “spreading wave”). Interestingly, the spreading wave of protests has indeed been observed in some of the events of social unrest [13,14].

There is a growing appreciation that the spatial aspect can shape the social protests in a number of ways; however, most of the previous modelling studies on protests’ dynamics focused on non-spatial or spatially implicit systems. Apparently, the consensus as to whether an adequate modelling framework should necessarily be spatial or can be non-spatial has not been reached yet. Whilst it is well known in other sciences that the dynamical properties of the spatially explicit system often are significantly different from those of its non-spatial counterpart [22], it remains unclear how the modelling predictions can differ in the context of the social dynamics. A more specific question is: How different can the dependence of the protester number on time be in spatial and non-spatial models for equivalent parameter values? Our paper intends to bridge this gap partially by considering a spatial extension of two non-spatial models recently developed in [17]. We will show that the dynamics of the spatial models can be qualitatively different, e.g., predicting the growth in the protester number for the parameter values where the non-spatial models predict decay. We will also show that the spatial model brings new dynamical features, in particular predicting a new scenario of the protests’ termination.

2. Mathematical Models in the Non-Spatial Case

We describe the dynamics of protests in terms of average values, hence disregarding any explicit effects of stochasticity: arguably, this approach works well when applied to sufficiently large systems [23], i.e., in the context of our study, to sufficiently large groups of people. We assume that the dynamics can be described by continuous, sufficiently smooth functions of time, so that in the non-spatial case, the mathematical model consists of ordinary differential equation(s).

We mention here that any group of people, e.g., the “population” of protesters or the larger community, is not uniform, but is structured according to various traits [15]. For instance, in the context of social unrest, the population of protesters can be structured according to the probability for a given individual to express a violent behaviour. Consider the Yellow Vest Movement as an example: whilst most participants protested peacefully, there was often a certain fraction that tended to behave aggressively. The wider community is structured with regard to the probability for a given individual to join the protests; this, in its turn, correlated with age, income, education, social background, etc. Whilst some people can be ignited easily (e.g., students [24,25]), people from other social groups might not so likely join the street actions in any circumstances. Aiming to build tractable models and to avoid unnecessary complexity (cf. “Occam’s razor”), in this paper, we disregard the effect of individual differences; in particular, we treat the protesters as a uniform group. We also assume that the total number of people in the society that could, under certain conditions, join the protests is sufficiently large and hence is not a factor that limits the growth of the number of protesters.

2.1. Single Species Model

We describe the civil unrest by a single variable, say N, which is the number of people actively participating in the street actions [26]. In order to build the model that describes the dynamics of N with time, we consider the following assumptions:

• The rate of change in the number of people attending the event is a result of the interplay between two processes, recruitment (people joining) and withdrawal (people quitting);
• Following earlier studies [18], we consider recruitment to be a collective phenomenon so that the recruitment rate is a nonlinear function of the number of people currently involved in the event;
• Decision of withdrawal is made individually, so that the withdrawal rate is a linear function of the number of people participating in the event, where the per capita withdrawal rate (say m) depends on time.

Under the above assumptions, we arrive at the following equation for the “population size” of the demonstrators:

$$\frac{dN\left(t\right)}{dt}=F\left(N\right)-m\left(t\right)N,$$

(1)

where the first and second terms in the right hand side are the recruitment rate and withdrawal rate, respectively, and $t\ge 0$ is time. We regard the per capita withdrawal rate m as a measure of the dissatisfaction with the effect of the protests, i.e., when the demands of the demonstrators are not being met by the authorities. In our model, higher dissatisfaction therefore means a higher rate at which people quit the movement. In case of an unchanging policy, the dissatisfaction can be expected to increase with time; hence, m is an increasing function of time.

We consider the recruitment rate in the following form:

$$F\left(N\right)={ϵ}_0+ϵN+\frac{a{N}^2}{h^2+N^2},$$

(2)

where ${ϵ}_0,ϵ,a$, and h are nonnegative parameters.

Parametrization (2) requires further discussion. We first notice that $F\left(0\right)={ϵ}_0\ge 0$. That reflects the observation that, in almost any society, there is often a fraction of people (albeit normally small) who are always ready and willing to rebel, e.g., to join a strike or a civil protest. Such people can be thought of as zealots that unconditionally support some ideas of ultimate social justice.

The per capita recruitment rate is given by the following expression:

$$\frac{1}{N}F\left(N\right)=\frac{{ϵ}_0}{N}+ϵ+\frac{aN}{h^2+N^2}.$$

(3)

For an intermediate range of N defined by conditions ${ϵ}_0/ϵ\ll N\ll h$ (if parameter values allow for the existence of such range), the per capita recruitment rate is therefore a monotonously increasing (approximately linear) function of N. This means that the rate at which new members are joining the movement tends to increase along with the number of people already in the movement. This is a typical behaviour response observed not only in humans [27], but also in many animals; in particular, it is believed to be one of the mechanisms resulting in the Allee effect [28]. However, when the number of people in the riots becomes large, $N\gg h$, the per capita recruitment rate is approximately constant, $\frac{1}{N}F\approx ϵ$. That reflects the assumption that there is a certain scale-independent probability of an average citizen to join protests, which agrees with observations [18,29].

Before considering the properties of Model (1) and (2) where the withdrawal rate m is a function of time, it is instructive to start with a simpler case where m is a constant parameter. In this case, it is readily seen that there exists two critical values, say $m_{\ast }$ and $m^{\ast }$, so that in an intermediate range $m_{\ast }<m<m^{\ast }$, the system possesses three steady states, two of them stable and the third one unstable; see Figure 1. The system undergoes a saddle-node bifurcation when m passes any of the two critical values; for $m<m_{\ast }$ and for $m>m^{\ast }$, there is only one steady state, large population state and small population state, respectively. In the general case ${ϵ}_0\ge 0$, an analytical expression for $m_{\ast }$ and $m^{\ast }$ is not available. In a special case ${ϵ}_0=0$, they are given as:

$$m_{\ast }=ϵ\mathrm{and}{m}^{\ast }=ϵ+\frac{\sqrt{2a-1}}{2h}.$$

In the case of a variable m, in order to derive the specific expression for $m\left(t\right)$, we take into account that, in the modern society, even in the absence of major events, a number of small scale protests are always happening. We assume that there is a certain background rate, say $m_0$, at which the participants of those protests quit the actions. We then further assume that the society is in a globally stable state, so that any perturbation from the background withdrawal rate should decay with time. Trying to keep the mathematical framework as simple as possible, the simplest equation that describes the corresponding evolution of m is given as follows:

$$\frac{dm\left(t\right)}{dt}=b(m_0-m),$$

(4)

where b is a coefficient. Equation (4) is readily solved, resulting in:

$$m\left(t\right)=m_1+(m_0-m_1)\left(1-e^{-bt}\right),$$

(5)

where $m_1=m\left(0\right)$ is the initial value of the withdrawal rate. In case the initial value is small, which can occur as a reflection of high social tensions, $m\left(t\right)$ is an increasing function. This seems to be in agreement with heuristic arguments. Note that, since m has the dimension of inverse time, it can be written as $m=1/\tau$ where $\tau$ is a characteristic duration of the individual participation. In a general case, $\tau$ can change in the course of the protests, hence being a function of time. At the beginning of the protests, the people’s expectations are high, the morale is high, and people are ready to spend a long time participating in the event. The characteristic duration of individual involvement is long, and the initial value of the withdrawal rate $m\left(0\right)=m_1$ is small. In the course of the protests, in case there is no sign of changes, the protesters’ morale may start falling, and the characteristic duration of individual participation could start decreasing.

Equations (1) with (2) and (5) are solved numerically. We consider $m_0$ and $m_1$ as the controlling parameters. Figure 2 shows the solutions obtained for different values of $m_0$ and $m_1$ and other parameters chosen as ${ϵ}_0=0.1$, $ϵ=3$, $a=146$, $h=35$, and $b=0.01$. (We mention here that our choice of parameter values was not entirely random or hypothetical: it was shown in [17] that, with this parameter set, Model (1)–(5) provided a good description of the Yellow Vests protests.) We therefore observe that there are solution properties that are rather robust to parameter variation. In particular, in all cases shown in Figure 2, the population size creates a peak (which can be rather flat for $m_1$ large enough; see Curves 5 and 6 in Figure 2b) followed by a gradual decay. This behaviour is in qualitative agreement with the dynamics of protests often observed in real life. The larger is $m_0$, the faster the rate of decay is. For a subcritical value of $m_0$, i.e., for $m_0<m^{\ast }$, the population size never falls to the low background value (which is zero for ${ϵ}_0=0$), but instead stabilizes at a higher positive value (cf. the magenta curve in Figure 2a), which corresponds to the upper stable steady state of the system; see the upper black circle in Figure 1b. For a supercritical value $m_0>m^{\ast }$, the upper steady state does not exist, so that in the large time limit N eventually converges to its lower steady state value (lower black circle in Figure 1b). However, if $0<m_0-m^{\ast }\ll 1$, the rate of convergence at intermediate time can be slow because of the “ghost attractor” [30,31]; e.g., see the middle part of the curves in Figure 2b.

2.2. Two Component Model

In the previous section, the withdrawal rate was described by the self-contained Equation (4) where the right hand side depends on m only. One weakness of this approach is that it completely disregards the possible dependence of the right hand side on N. Even in a more general case where the right hand side in Equation (4) could be a nonlinear function, m is effectively just a given function of time. As an immediate yet nontrivial extension of the above model, we now consider the case where the population dynamics of riots is described by the following equations:

$$\frac{dN\left(t\right)}{dt}=F\left(N\right)-mN,\frac{dm\left(t\right)}{dt}=G(m,N),$$

(6)

where G is a certain function to be specified.

Apparently, Model (6) is more general than Model (1) and (2) and hence can describe a broader variety of situations. In order to find an appropriate parametrization of G, we have to consider more closely the reasons why people may quit the protests. In doing this, we first observe that street actions usually bring some economic damage. That can include direct damage (such as broken shop windows, burnt cars, etc.) and indirect damage to the economy (e.g., through the disruption of public transport). We assume that the amount of damage caused is proportionate to the number of people involved in the street actions. We then further assume that an ordinary protest participant eventually develops the feeling of guilt for causing damage to the community and that this feeling is likely to cut short the duration of his/her participation in the street actions. Correspondingly, this means that the rate of change in m should be an increasing function of the damage, which we consider to be linear.

Correspondingly, we now arrive at the following model of the protests’ dynamics:

$$\begin{array}{c}\hfill \frac{dN}{dt}=\left({ϵ}_0+ϵN+\frac{a{N}^2}{h^2+N^2}\right)-mN,\frac{dm}{dt}=\beta N,\end{array}$$

(7)

where $\beta$ is a coefficient. Unless the protests are extremely violent, it seems reasonable to assume that $\beta \ll 1$.

Figure 3 shows the phase plane of Model (7). It is readily seen that the properties of Model (7) are significantly different from those of Model (1) and (2). In particular, there is no steady state in the interim of the first quarter of the plane, i.e., for $N>0$ and $m>0$.

In order to reveal the solution properties, the Equations (7) are solved numerically. We consider $\beta$ as the controlling parameter. Typical results are shown in Figure 4. At the beginning, the number of protesters promptly becomes very large, thus creating a high peak (see Figure 4a). Past the peak, the number of protesters starts falling, first fast, then at a decreasing rate. The number of protesters keeps decreasing with time at a lower and lower rate until, at a certain moment, it promptly decays to zero. The corresponding dependence $m\left(t\right)$ is shown in Figure 4b.

The behaviour of the solution of Model (7) is thus apparently similar to that of the solution of the single component Model (1) and (2); cf. Figure 2. There are, however, some important differences. In Model (1) and (2), there exists a finite critical value $m^{\ast }$, so that for $m_0<m^{\ast }$, the number of protesters never falls to the small equilibrium value, but stabilizes at the higher steady state (see Curve 1 in Figure 2a). On the contrary, in Model (7), there is no finite critical value: for any $\beta >0$, the solution $N\left(t\right)$ will eventually converge to the lower steady state. We also note here that the shape of the graph of $N\left(t\right)$ shown in Figure 4a, albeit similar to Figure 2, is not the effect of the ghost attractor, but a property of the so-called slow fast dynamics (when $\beta \ll 1$) [31].

3. Spatially Explicit Model

In the previous section, we completely disregarded the existence of the spatial dimension in the dynamics of the social protests, hence implicitly assuming that the population density of the protesters is distributed uniformly over space. In reality, it is rarely so, especially when the dynamics is considered over a sufficiently large spatial scale such as a region, a county, a province, or the whole state. There are several reasons for the spatial heterogeneity to emerge. Firstly, the overall population density itself is distributed highly non-uniformly, i.e., reaching very high values in urban areas, but falling to low values in the countryside. Secondly, people’s behavioural response and their social activity can be significantly different at different locations in space, e.g., as a result of different living conditions and different accessibility to social media. Therefore, although the consideration of non-spatial models makes the first necessary step in the study, a more adequate mathematical framework for modelling social protests must take space into account explicitly, as has indeed been recognized in a number of studies, e.g., see [8,12,13,14]. In this section, we consider the models of the social protests dynamics in space and time that consist of reaction-diffusion equations, where the “reaction” part is given by the corresponding non-spatial model (as in Section 2) and the spread of the protests in space (described by the diffusion terms) is regarded as a certain “contagion” [13,32].

We mention here that the spatial spread of information and emotions and hence the spread of social tensions are known to have two different spatial modes [12,32]. In one mode, people communicate in direct encounters due to their geographical proximity: arguably, people living and/or working in the same area would normally meet in person more frequently than people living far apart. The importance of this mode of contagion is indirectly confirmed by the well established fact that most of everyday human travel occurs within a short range from peoples’ homes [33]. In the second mode, the spread of information and emotions is not related to travel and occurs due to the connection by social media. For the first mode, the spread of protests can be described by the standard diffusion [12,14]; for the second mode, kernel based integro-differential equations are thought to be more adequate [13] as they better describe the effect of long distance connections [34]; see also [35,36,37].

Note that, whilst the effect of the first mode is obviously limited to short distances, the second mode is not restricted to any particular spatial scale: indeed, technically, social media can be accessed from anywhere in the world. Hence, one can think, respectively, about the local coupling vs. the global coupling. In the case of the local coupling due to direct encounters, the effect of space is likely to be more pronounced; in particular, the level of social tensions should be a function of space. On the contrary, in the case of the global coupling due to social media, the social tension is likely to be distributed over space more uniformly.

3.1. Single Species Model

We begin with the single species system. In terms of the above discussion, it corresponds to the global coupling where the social tension is distributed over space uniformly, and hence, the withdrawal rate does not depend on space (but depends on time). The distribution of the “population” of protesters over space is then described by a single dynamical variable, i.e., by the population density, say $u(x,t)$. The spatially explicit counterpart of Model (1) and (2) is described, in the 1D case, by the following equation:

$$\frac{\partial u(x,t)}{\partial t}=\left({ϵ}_0+ϵu+\frac{a{u}^2}{h^2+u^2}\right)-m\left(t\right)u+D\frac{{\partial }^{2}u}{\partial x^2},$$

(8)

where x is the position in space, $-L\le x\le L$, D is the diffusion coefficient, and the dependence of the withdrawal rate m on time is given by Equation (5).

We assume that protests start at a certain location or a small area. Correspondingly, for the initial conditions, we consider a piece-wise constant distribution of finite support described as follows:

$$\begin{array}{c}\hfill u(x,0)=u_0\mathrm{for}-\Delta \le x\le \Delta ,u(x,0)=0\mathrm{otherwise}.\end{array}$$

(9)

At the ends of the domain, we use the Dirichlet boundary condition, $u(-L,t)=u(L,t)=0$.

Equation (8) with (9) was solved numerically by finite differences. Typical results are shown in Figure 5. We readily observe that the system dynamics are characteristic for reaction-diffusion systems with compact initial conditions, albeit being modified by the dependence of the withdrawal rate on time. At an early stage of the dynamics, the evolution of the initial distribution results in the formation of two travelling fronts propagating away from the location of the initial condition (see the yellow, violet, green, and light-blue curves in Figure 5a): the protests spread fast over space. In the course of time, the speed of the fronts decreases as a result of the increase in m. Eventually, the front stops and reverses (cf. the dark-blue and brown curves in Figure 5a). The population then retreats to the place of its original introduction, which finally results in its collapse and extinction, i.e., the termination of the protests.

We mention here that we are not aware of any well developed mathematical theory of travelling waves in non-autonomous reaction-diffusion systems, i.e., with reaction rates explicitly depending on time. It seems, however, possible to explain the numerical results in terms of autonomous systems. In the baseline case of constant m, the reaction term of Equation (8), that is:

$$\begin{array}{c}\hfill \mathcal{F}\left(u\right)={ϵ}_0+ϵu+\frac{a{u}^2}{h^2+u^2}-mu,\end{array}$$

(10)

is a sigmoidal curve, which, in the general case $m_{\ast }<m<m^{\ast }$ (cf. Figure 1), has three roots, $0\le u_{lower}<u_{int}<u_{higher}$ (where $u_{lower}=0$ if ${ϵ}_0=0$). The direction of the front propagation is known to depend on the value of the following integral [38,39]:

$$\begin{array}{c}\hfill M={\int }_{u_{lower}}^{u_{higher}}\mathcal{F}\left(u\right)du,\end{array}$$

(11)

so that the front propagates towards the low density area (and hence away from the place of the compact initial condition) if $M>0$ and towards the high density area if $M<0$. In the former case, the protests invade, and in the latter case, the protests retreat.

In our case, m depends on time, and strictly speaking, the above results are not applicable. However, one can expect that they may still work approximately. Indeed, it is readily seen that, for the parameters chosen as in Figure 5, $M>0$ at an early time, but becomes negative at a large time, changing its sign at $t\approx 123$ (see Figure 6). This explains the reversing of the front.

We therefore conclude that the spatial dynamics exhibits some new properties that did not exist in its non-spatial counterpart, in particular by providing a new scenario of the termination of the social unrest. Recall that, in the non-spatial system, the unrest terminates only after the increase in m results in the disappearance of the upper positive steady state (see Figure 1). In the spatial system, the unrest termination can happen even when the positive steady state exists: for the termination to occur, it is sufficient that m becomes sufficiently large (but not necessarily supercritical) to reverse the front propagation.

The dependence of the total number of protesters on time, i.e., $N_{tot}\left(t\right)={\int }_{[-L,L]}u(x,t)dx$, appears to be rather different in the spatial system compared to the corresponding non-spatial system; see Figure 7. Instead of the sharp peak, after the fast increase at the beginning, $N\left(t\right)$ then exhibits a much slower decrease than in the non-spatial system (cf. Figure 2), eventually accelerating to a fast decay to zero.

3.2. Two Component Model

We now consider the spatially explicit version of Model (7), which is given by the following equations:

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& =& {ϵ}_0+(ϵ-m)u+\frac{a{u}^2}{h^2+u^2}+D\frac{{\partial }^{2}u}{\partial x^2},\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \frac{\partial m(x,t)}{\partial t}& =& \beta u+D_m\frac{{\partial }^{2}m}{\partial x^2}.\hfill \end{array}$$

(13)

For the sake of generality, we assume that the spatial dynamics of the withdrawal rate may occur with a different diffusion coefficient $D_m$.

For the initial conditions, similarly to the above, we consider the piece-wise constant distributions:

$$\begin{array}{cc}\hfill & u(x,0)=u_0\mathrm{for}-\Delta \le x\le \Delta ,\mathrm{otherwise}u(x,0)=0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & m(x,0)=m_0\mathrm{for}-{\Delta }_m\le x\le {\Delta }_m,\mathrm{otherwise}m(x,0)=0,\hfill \end{array}$$

(15)

where ${\Delta }_m\le \Delta$ to reflect the assumption that there is no damage without protests. At the ends of the domain, we use the Dirichlet boundary condition, $u(-L,t)=u(L,t)=0$ and $m(-L,t)=m(L,t)=0$.

Mathematical Problem (12)–(15) was solved numerically by finite differences using the forward scheme. The mesh steps were chosen sufficiently small to avoid numerical artefacts, and it was checked that a further decrease of the steps value had no effect on the numerical solution.

Some typical results are shown in Figure 8. We readily observe that, as well as in the previous case, the evolution of the initial conditions results in the formation of two travelling waves propagating away from the location of the initial distribution. However, the properties of the wave are essentially different from those in the single species model. Instead of a smooth, monotonous front, we now have a travelling peak with a very high population density at its maximum. In the wake of the peak, the population density promptly decays to a much smaller value. The corresponding spatial distribution of the withdrawal rate is shown in Figure 9. After the peak reaches the domain boundary, it promptly disappears, and the total number of the protesters decays to zero.

As we observed in numerous numerical simulations, the main properties of the solution, such as the formation of the population peaks propagating away from the location of the initial distribution, were generic and robust to parameter changes. In particular, we did not observe any significant changes in the system dynamics for different values of the diffusivity ratio $D_m/D$ (having checked several values of the ratio in the range between one and 0.01). However, some quantitative features of the dynamics appeared to be somewhat sensitive to other parameter values. For instance, the height of the peak depended on $ϵ$ rather strongly. An example shown in Figure 10 is obtained for $ϵ=0.1$. We readily observe that, whilst the dynamics is qualitatively the same as in the case $ϵ=1$ (see Figure 8), the height of the peak is almost an order of magnitude smaller. The speed of the peak propagation is slightly less as well. The distribution of damage, however, shows features similar to what was obtained for $ϵ=1$ (see Figure 9); we do not show it here for the sake of brevity.

Figure 11 shows the total number of protesters vs. time as is given by $N_{tot}\left(t\right)={\int }_{[-L,L]}u(x,t)dx$. It is readily seen that the dynamics of the protests in the spatially explicit Model (12) and (13) is considerably different from its non-spatial counterpart (7). Recall that in the non-spatial system, the number of protesters reaches its maximum at early time and then decays monotonously, first very fast and later at a decreasing rate (see Figure 4a). In the spatial system, however, the fast increase at the beginning changes to a slower growth that goes on for a considerable length of time. The phase of slow growth corresponds to the peak travelling over space where the shape of the peak remains approximately the same, but the length of its tail grows steadily. The duration of this phase depends on the extent of the available space; the larger the spatial domain is, the longer the total number of protesters keeps growing slowly, but steadily. After the peak hits the boundary of the domain, the effect of the boundary conditions bring it down, resulting in the sharp decrease in the protesters’ number.

Note that the collapse of the travelling wave as a result of the collision with the boundary is a common phenomenon in reaction-diffusion systems of various origins [39,40,41]. Such a collapse is not very sensitive to the type of the boundary condition; in particular, in our model, it was observed both for the Dirichlet and Neumann types (for the sake of brevity, we do not show here the corresponding figures for the zero-flux boundary conditions). We also mention that the collapse of protests observed in our model when the wave reached the domain boundary could be insightful for understanding the dynamics of real events of social unrest. For example, in the 2005 French riots, it was reported that the street protests were quickly spreading across suburbs before they reached some quieter areas, after which the whole unrest went extinct [13]. Although the collapse of the 2005 unrest in Paris can be explained by a number of reasons, e.g., overall decrease of the protesters’ motivation, increasing efficiency of the police, etc., reaching the edge of the area where the social or geographical conditions changes significantly is arguably a major factor.

Apparently, the rate of increase in the total number of protesters depends on the speed of the travelling peak. That raises an interesting mathematical question as to whether the speed can be estimated. We recall here that, in reaction-diffusion systems, the travelling wave usually arises as a heteroclinic connection between two equilibria (homoclinic in the case of a travelling peak) in the phase space of the system, the speed of the travelling wave being an eigenvalue of the corresponding nonlinear problem [39,42,43]. Our system (12) and (13) does not quite fit into this framework, as a positive equilibrium does not exist for ${ϵ}_0=0$. This indicates that, in a strict mathematical sense, the travelling wave does not exist. However, arguably, the part of the solution that describes the peak still behaves as a travelling wave, as indeed is readily seen from the properties of the numerical solution, e.g., see Figure 8 and Figure 10. Hence, we hypothesize that some of the tools normally used to estimate the speed of the wave may still be applicable.

Such tools must be applied with due care though. We first observe that the solution behaves as a travelling wave at the “leading edge” of the travelling peak, i.e., where both components u and m are very small. Equation (12) can therefore be linearized taking the following very simple form:

$$\begin{array}{c}\hfill \frac{\partial u(x,t)}{\partial t}=ϵu+D\frac{{\partial }^{2}u}{\partial x^2},\end{array}$$

(16)

where we consider ${ϵ}_0=0$ here and below for the sake of simplicity. Note that the linearized equation does not contain m any more. For the linear reaction-diffusion equation with finite or semi-finite initial condition, it is well known that, although the travelling wave solution does not exist globally, the leading edge of the solution propagates as a travelling wave with the speed given by:

$$\begin{array}{c}\hfill c=2\sqrt{ϵD},\end{array}$$

(17)

e.g., see [37,38]. Estimate (17) obtained by linearizing the system (12) and (13) therefore coincides with the speed of the travelling front in a scalar reaction-diffusion equation with the generalized logistic growth [42,43].

However, the expression (17) appears to be in bad disagreement with numerical results predicting a value considerably less than was obtained in the simulations, e.g., two instead of 3.1 for the parameters of Figure 8. Even worse, it does not in any way predict the dependence of the speed on parameter values: for instance, for $ϵ=0.1$ Expression (17) gives $c\approx 0.63$, whilst the actual value is about 2.6. A more careful look readily reveals the reason for this disagreement: as a matter of fact, Expression (17) is hardly relevant at all as its validity is limited to the case where the growth function is convex [42,43], whilst in our case, it is not (see Figure 12).

We therefore have to apply a different approach, namely the approach based on the comparison principle for PDEs [44]. We first notice that Equation (12) can be written as the following differential inequality:

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& \le & ϵu+\frac{a{u}^2}{h^2+u^2}+D\frac{{\partial }^{2}u}{\partial x^2}.\hfill \end{array}$$

(18)

We now build an auxiliary function $\tilde{\mathcal{F}}\left(u\right)$ as follows:

$$\begin{array}{ccc}\hfill \tilde{\mathcal{F}}\left(u\right)& =& \alpha u\mathrm{for}0\le u\le u_A,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill \tilde{\mathcal{F}}\left(u\right)& =& ϵu+\frac{a{u}^2}{h^2+u^2}\mathrm{for}u>u_A,\hfill \end{array}$$

(20)

where $\alpha$ is the slope of the tangent line to the graph of the reaction term and $u_A$ is the abscissa of the conjunction point (Figure 12). After some simple algebra, we obtain that:

$$\begin{array}{c}\hfill \alpha =\frac{a}{2h}+ϵ.\end{array}$$

(21)

Obviously, $\tilde{\mathcal{F}}\left(u\right)$ is an upper bound for the reaction term in the right hand side of (18), so that the following inequality holds:

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& \le & \tilde{\mathcal{F}}\left(u\right)+D\frac{{\partial }^{2}u}{\partial x^2}.\hfill \end{array}$$

(22)

By the virtue of the comparison principle, the solution of the following equation:

$$\begin{array}{ccc}\hfill \frac{\partial \tilde{u}(x,t)}{\partial t}& =& \tilde{\mathcal{F}}\left(\tilde{u}\right)+D\frac{{\partial }^2\tilde{u}}{\partial x^2}.\hfill \end{array}$$

(23)

gives an upper bound for the solution $u(x,t)$ of Equation (12) provided the initial and boundary conditions are the same.

Function $\tilde{\mathcal{F}}\left(u\right)$ is convex; therefore, the linearization at the leading edge can be applied resulting, taking (21) into account, in the following upper bound for the speed of the wave:

$$\begin{array}{c}\hfill c=2\sqrt{\alpha D}=2\sqrt{\left(\frac{a}{2h}+ϵ\right)D}.\end{array}$$

(24)

Figure 13 shows the values of the speed calculated in simulations along with the analytical estimate (24). It is readily seen that, while the expression (24) somewhat overestimates the actual value, it manages to describes reasonably well the dependence of the speed on the parameter, i.e., the approximately linear increase along with an increase in $ϵ$.

4. Discussion and Concluding Remarks

Social protests are a complex phenomenon that has a variety of important implications for the functioning of different social groups and organizations, as well as the society as a whole. Hence, a good understanding of the protests’ dynamics is needed, in particular to identify the factors controlling their intensity and duration and to reveal the generic patterns of the protests’ evolution in time and space. Whilst this is apparently a problem in the remit of social sciences, mathematical modelling has long been regarded as a powerful research tool that helps to achieve these goals [3,18].

In this paper, we considered a few generic, conceptual models of the social protests dynamics. We disregarded the social structure of the crowd and described the dynamics by a single dynamic variable, i.e., the number of protesters. Our study was partially motivated by the observation that there can be two apparently different scenarios of the protests’ development in time: the one where the peak in the total number of protesters is followed by a fast decay and the other one where the decay in the wake of the main peak is very slow, thus creating a long tail. Correspondingly, our first goal was to understand whether these two apparently different cases can be embraced in the framework of the same model. We found that the two cases differed quantitatively rather than qualitatively, i.e., the different decay rate (and hence, the total duration of the protests) was determined by the model parameters. Interestingly, we obtained that the duration of the protests could, in principle, be indefinitely long, although the corresponding parameter range was small.

There have been several studies previously where riots and social unrest were treated in terms of nonlinear dynamical systems [2,11,12,13,17]. Adopting this modelling paradigm, the slow decay in the protests intensity could be linked to the so-called long transient dynamics (cf. [17]): a phenomenon inherent in many dynamical systems of different origin [30,31,45,46]. In a nutshell, a long transient means that the dynamical system mimics its stable asymptotical dynamics (e.g., steady state, periodical, or chaotic), being in fact in a transient regime. The quasi-asymptotical dynamics can persist for an indefinitely long time before the system experiences a fast transition to its true asymptotics. In a non-spatial system, the main mechanisms resulting in long transient are ghost attractors, crawl-bys, and slow-fast dynamics [30,31], and this is precisely what determines that long tail predicted by our models. The long tail obtained in the single species model (cf. Curves 1–3 in Figure 2b) was the effect of the ghost attractor that emerged after the upper stable steady state of the system disappeared in a saddle-node bifurcation when the final value of the withdrawal rate exceeded the critical value $m^{\ast }$ (see Figure 1b). The long tail predicted by the two component model (e.g., see the green curve in Figure 4a) was a manifestation of the slow-fast dynamics that arose for small values of $\beta$.

We mention here that linking our findings to some general properties of the dynamical systems (such as the existence of long transient dynamics) is not a futile mathematical exercise, but brings insights that can have immediate practical value. In particular, the generic property of long transients that the slow decay in the protesters number is followed by a fast decay to zero may be used as a basis for a response strategy: if the authorities have enough patience to wait (or other reasons not to use excessive force in response to protests), it is likely that the protests eventually disappear due to their inherent dynamics (perhaps this is what is happening with the 2019 riots in Hong Kong where the authorities have been surprisingly mild in their response, probably just waiting until the riots burn out by themselves).

The second goal of this study was to highlight the role of the spatial dimension in the dynamics of social protests. We used mathematical modelling to reveal essential differences between the properties of the spatial systems and their non-spatial counterparts. The importance of space is well recognized in natural systems, in particular in ecology [22], but its effect on the social processes remains somewhat obscure. We obtained that both models predicted the emergence of a travelling wave of the protests’ activity (albeit with rather different properties). This was in agreement with other recent modelling studies [13,14]; also, the spatial dynamics reminiscent of a spreading wave was observed in some data [13,47].

Our results suggested two essentially different scenarios of the protests’ development in space, say, Type I and Type II. In the Type I scenario, the protesters are distributed over a certain area with approximately uniform density, this area initially growing with time. The growth of the area “invaded” by the protests is combined with a slow decrease in time (apart from a very early stage of the dynamics) of the protesters’ population density uniformly over the whole invaded space, eventually resulting in the protests shrinking back to the place of their origin and collapsing (see Figure 5). This is a result of the global coupling where the withdrawal rate is the same everywhere in space: m depends on time, but not on space. In the Type II scenario, the distribution of protesters over space is highly heterogeneous. The large density of protesters is only observed inside a narrow range or small area, which in the course of time moves away from the place where the protests started (see Figure 8). The collapse of the protests occurs as a result of a specific “boundary forcing”, i.e., when the peaks are pushed close to the domain boundary. This is a result of the interplay between the local coupling (where, generally speaking, the withdrawal rate is different at different locations) and the feedback mechanism through which the withdrawal rate depends on the protesters’ density. One can further speculate as to how the difference between the two scenarios (and between the models, respectively) can be projected onto real life. Since the explicit dependence of the withdrawal rate on space may indicate the heterogeneity in the peoples’ responses, one can suggest that the Type I scenario is more likely to occur in a well connected (hence global coupling), socially uniform society, while the Type II scenario is more likely to be seen in a poorly connected society with a broader range of social norms.

Note that both non-spatial models that we considered in Section 2 describe a very similar pattern in time, i.e., the main peak was followed by a tail, the rate of decay at the tail depending on the model parameters. Since different models were built to account for the effect of different factors, that evokes the old problem of the relation between the pattern and the process: if two different models predict the same pattern, how do we “let the right one in”, i.e., what are the processes that are actually important? However, when the models are extended to include space, the dynamics become qualitatively different. Therefore, here, we showed that taking into account the spatial dimension of the system could be one possible way to solve this problem.

Our study leaves a few open questions. Firstly, in our models, we assumed that the “population” of protesters consisted of essentially identical individuals. In reality, people are different in terms of beliefs, responses, social status, etc. [15,18,27,48,49]. How the heterogeneity of the population can modify our findings remains unclear. In particular, the recruitment of protesters can be facilitated by activists, as was the case in the Russian revolutions in 1905–1907 and later in 1917 [1,19,20]. Currently, such activists are active Internet users and social engineers and may have even a stronger influence. Secondly, with regard to the spatial extension of the models, we restricted our analysis to the hypothetical 1D case. How the patterns of spread may change in the more realistic 2D system remains to be seen. Heterogeneity of space can also be an important factor affecting the dynamics of protests, as was shown in a number of works [9,12,13], and that may require a mathematical framework different from the standard Fickian diffusion, e.g., networks. Consideration of these factors should become a focus of future research.