ANALYSIS AND SIMULATIONS OF THE HANDY MODEL WITH SOCIAL MOBILITY, RENEWABLES AND NONRENEWABLES

. We expand the HANDY (Human And Nature DYnamics) model for the socioeconomic dynamics of a large stratiﬁed society. The basic model was introduced in Motesharrei (Dissertation 2014) and Motesharrei et al. (2016). It is a nonlinear system of ODEs for a ‘simple society’ of Elites, Workers, Wealth, and Natural Resources. Following Ali Kadhim (Dissertation 2021), we add social mobility between the classes and split natural resources into renewables and nonrenewables. We establish the existence, boundedness and positivity of the solutions, and investigate the stability of the steady states. The model admits stable steady states, and there is numerical evidence of stable periodic solutions and limit cycles. Simulations depict the diﬀerent qualitative types of model behavior: convergence to steady states, periodic oscillations or collapse.


Introduction
The HANDY -Human And Nature DYnamics -was a ground breaking sustainability focused mathematical model for complex societies, constructed in Motesharrei's dissertation [5].The original model consisted of four ordinary differential equations (ODEs), representing the growth rates of two populations, their use of natural resources and of their wealth.The model was simulated to study its behavior in [5] and in Motesharrei et al. [6,7].Computer simulations of the HANDY model showed that economic stratification and rapid depletion of natural resources were among the main reasons for societal collapse.The simulation results indicated that an irreversible collapse could be avoided if the society adopted sustainable behaviors to reduce the depletion of natural resources and moved towards an equitable distribution of resources.In these publications, the authors provided many historical examples and studies of collapse of earlier civilizations, looking for possible underlying processes that were common to all, instead of trying to study each society separately.Indeed, the main interest in this 'simple model' and its computer simulations was to find common causes for the collapse of the civilizations, in addition to issues specific to each one.It was found that the common causes in all the cases were both over-depletion of resources and inequality.This finding of common causes for the collapse of the different societies is what makes mathematical models so important and effective in clarifying some general types of behavior that are hard to see when studying each society separately.
Related previous work is as follows.In Brendan and Taylor [2], the authors introduced a simpler model to study the collapse of the society in Easter Island.They used the Ricardo-Malthus model with two ODE's, one for renewable resources that included a logistic term with saturation and a harvesting term.The population growth rate was described by three terms: the birth rate, the death rate, and a fertility factor using the Malthusian term.The latter claims that higher income rate stimulates population growth, leading to the depletion of resources and eventually to famine.On the other hand, it was found that too low wages may cause a population collapse caused by conflict.Various civilizations have gone through such collapse, including the western Roman empire and the Mayan and the Mesopotamian civilizations.A similar model for the societal dynamics was considered by Grammaticos et al. in [3], consisting of three ODEs representing the population, resources and reserves.They used a special factor for the amount of reserve per capita, and proposed that the population was flourishing when the accumulated reserve was abundant (above the threshold) and died out when the reserves were insufficient.Furthermore, they found in their simulations patterns, similar to those in the HANDY model, of nature's slow growth and rapid decline, known as the Seneca effect.
This work continues the study of such models, and is based on the dissertation [1], where a few different versions of the model were constructed, analyzed and computer simulated.Here, we have a three-fold aim: First, we expand the HANDY model by allowing for social mobility between the workers and the elite, and by splitting the resources into renewable energy/wealth sources (wood, solar, wind), and nonrenewable sources (oil, gas, coal).Secondly, we perform mathematical analysis of the model, proving the existence, boundedness and positivity of the solutions, as well as a study of the system steady states, following [1].Thirdly, we present simulations of typical types of model predictions, and initial study of the dependence of the solutions on the consumption rates of the resources and on inequity in consumption.
The computer simulations indicate that the system has four qualitatively different types of behavior: monotone approach to the steady state, Figure 1 (a) (where the approach is monotone after one initial oscillation); damped oscillations approach, Figure 1 (b), periodic oscillations, Figure 1 (c), and system collapse, Figure 1 (d).The different cases depend on the different sets of initial conditions and the values of the various system parameters, and provide insight into possible behaviors of such models and hence, to the potential behavior of complex societies.It seems that such extensions of the HANDY model may be useful in investigating various scenarios and the associated conditions for flourishing or collapse of complex societies.It is very likely that the system is rich enough to exhibit chaotic solutions, however, our simulations did not find any.Moreover, we do not have a mathematical proof of the existence of periodic solutions, yet.
Section 3 establishes the existence, positivity and boundedness of the solutions of the model.Similar properties of the HANDY model have been established in [1], and the analysis here is tailored to the HANDY-SMRN model.Section 4, studies the system's steady states and their stability.It is found that the steady states depend essentially on the ratio η = w/w th of the wealth w to the wealth's threshold w th .Because of the structure of the system coefficients, there are three regimes, which lead to three families of steady states: Case (i), 0 ≤ η ≤ 1/κ 1 , which includes the two steady states with collapse, namely (0, 0, 0, 0, 0, ) and (0, 0, λ, 0, 0).The first of which is found to be unstable while the second is stable and attracting.Case (ii), 1/κ 1 ≤ η < 1, where the elite have enough wealth; and Case (iii), 1 < η, where there is enough wealth (which includes food reserves) for both populations, to avoid collapse.
We note that the system has 15 parameters, and five initial conditions, so a complete numerical study is out of reach.Therefore, we present only the four qualitatively different types of solutions; simulations with different δ n and δ r ; and with different κ 1 .We note that a more thorough computer study of the effects of different rates of depletion of renewables and nonrenewables may be of interest.
The computer simulations with different values of the depletion of nonrenewables, with two different values of δ n show that a more moderate use of nonrenewables leads to longer periods of reaching the steady states but also to substantially higher steady state populations.Concerning the renewable natural resources, with depletion factor δ r , we found that there is a range of values for which the solutions are either oscillatory, or oscillate for very long periods of time and then quickly converge to the steady states.
Finally, simulations with different consumption rate constant κ 1 show that for κ 1 = 1, 5, 10 the approach to the steady states is increasingly more complex with more large oscillations, in which the populations almost die out.
Note: It is an important observation that even in the cases where, eventually, the solutions converge to similar steady states, the ways these are reached may be very different, with more or less societal upheaval on the way.This indicates that just studying the steady states of the model and their stability may provide only a very partial and limited description of the society's dynamics.
Clearly, there is considerable interest to find out if the model can fit some of the historical collapses of societies which essentially consisted of commoners and nobility.
The rest of the work is organized as follows.The model is presented in Section 2, where the underlying assumptions are discussed with some detail, taken from the original HANDY model and [1].Then, we discuss our expansion of the model by allowing for mobility between the two classes, and splitting the natural resources into renewables and nonrenewables.We also introduce the ratio of wealth to the threshold wealth, η, which plays an important role in the analysis.Section 3 deals with the analysis of the model.It establishes the boundedness and positivity of the solutions, Proposition 3.2.These lead to the existence of the solutions on every finite time interval, Theorem 3.3.Then, Section 4 studies the system steady states and establishes the conditions for their stability.Section 5 presents the results of model simulations.It depict the four main types of behavior that the model predicts.Also it presents a short study of the effects of increasing or decreasing the renewables and nonrenewables rate constants, and the dependence on the consumption rate constant.Section 6 poses some possible unresolved questions of interest.

HANDY-SMRN model
The extended model studied here, is a version of the model constructed in Ali Kadhim [1] (Dissertation 2021; Chs. 5 and 6) and consists of five coupled nonlinear ODEs: two equations model the populations growth rates of the Commoners or workers, x c (t), and the Elites or the rich, x e (t); one describes the rate of growth or depletion of renewable natural resources, y r (t) (such as wood, hydro-power, or solar energy), and one for the depletion of non-renewables (such as coal, oil or gas) y n (t), and the rate of growth of wealth w(t), which includes food surpluses.The time t is measured in 'time units,' (in real societies this would be in years).
The novelties in this model, as noted in the introduction, are the addition of social mobility between the rich and the workers, and the splitting of the natural resources into renewables and nonrenewables.
The complete HANDY-SMRN model, the HANDY model with Social Mobility, Renewables and Nonrenewables, is the following: Problem 2.1.Find five functions (x c (t), x e (t), y r (t), y n (t), w(t)), defined on [0, T ], 0 < T < ∞, such that ) ) Here, the consumption rate functions C c and C e are given in (2.9), the wealth threshold w th is given in (2.7), and the death rates α c and α e are given in (2.10).The initial conditions are, x c (0) = x c0 , x e (0) = x e0 , y r (0) = y r0 , y n (0) = y n0 w(0) = w 0 . (2.6) We assume that x c0 > 0, x e0 ≥ 0, y r0 > 0, y n0 ≥ 0, w 0 ≥ 0, are prescribed.The rest of the notation is as follows.The respective birth rates are β c , β e , assumed to be constant; γ is the nature's regeneration factor, λ its saturation level or the carrying capacity of the renewable resources; δ r and δ n are the respective renewable and non-renewable resources depletion rate constants; δ w is the conversion factor between the resources depletion and wealth generation.Next, γ e is the rate constant in which the rich become workers or commoners (going bankrupt), and γ c is the rate constant in which the workers 'make it' and become wealthy.All these may be constants or given non-negative functions.
We note that (2.3) is of the logistic type, and allows for regeneration of the renewable resources up to the carrying capacity λ, while equation (2.4) for the nonrenewable does not allow for regeneration.
Before we describe C c , C e , the total consumption rates, we introduce the wealth threshold where > 0 is a very small number that prevents the threshold w th from becoming zero, since it appears in the denominators in (2.9) and (2.10).We next introduce the variable which plays an important role in the analysis of the steady states, (2.8) The total consumption rates are given by (2.9) Finally, the death rates α c , α e , are given by, , (2.10) We refer to [5] and to [1] for a more detailed description of the model parameters and coefficients.We note that in the original HANDY model [5], the choice was

Model analysis
This section provides the analysis of the model, starting with establishing the boundedness, positivity and then the existence of the solutions to Problem 2.1.Then, we study the steady states of the system and establish their stability.
The existence of the solution is based on the following theorem that can be found in Kuttler [4], and for the benefit of the reader, its statement is included.Theorem 3.1.Let F : [0, T ] × R n → R n be continuous and suppose that there exists L > 0 such that for all λ ∈ (0, 1), if for all t ∈ [0, T ], implies that x < L.Then, there exists a solution to the system of ODEs, We note that the theorem does not guarantee the uniqueness of the solution.To show uniqueness, a more refined argument is needed, which we leave open.

3.1.
Positivity and a priori estimates.Let z = (x c , x e , y r , y n , w) and F : R 5 → R 5 be given by Here, C c and C e are given in (2.9), w th is given in (2.7), and α c and α e are given in (2.10).We assume, following [5,1], that 0 < α m < β e ≤ β c ≤ α M < 1, and then it follows that Everywhere below, we assume that κ 0 , κ 1 ≥ 1.
First, we write the system as Here and below, the prime indicates the time derivative, and each component of z 0 is assumed to be positive.Below, we show that the results hold true for x e0 = 0, y n = 0 and w 0 = 0, as well.
We next establish the necessary a-priori estimates on the system solutions.Let 0 < T < ∞, be fixed, and then [0, T ] is an arbitrary finite interval.We begin by showing that starting with positive initial conditions, all the possible solutions are nonnegative and bounded.This leads to estimates that are valid for [0, T ] and allow us to prove that a solution exists on this interval.Proposition 3.2.Assume that x c0 > 0, y r0 > 0, y n0 > 0 and x e0 , w 0 are nonnegative.If z = (x c , x e , y r , y n , w) is any solution to the system (3.5), then all the components of z are nonnegative.Moreover, for 0 ≤ t ≤ T , the following estimates hold: Here, M is a positive constant (see below) that depends only on the problem parameters.
Proof.It follows from the first two equations, by multiplying the first with x e and the second with x c , adding the expressions and rearranging, that e .Since initially x e (0)x c (0) > 0, by continuity, there is t 0 > 0 such that x e (t)x c (t) > 0 for 0 ≤ t < t 0 .Then, if x e (t 0 )x c (t 0 ) = 0, we obtain e ≥ 0. If either x e (t 0 ) > 0 or x c (t 0 ) > 0, then (x e x c ) (t 0 ) > 0, so, by continuity, the function (x e x c )(t) is increasing at t 0 , which contradicts the assumption that x e (t 0 )x c (t 0 ) = 0 for the first time.Now, suppose that both vanish at t 0 > 0.Then, for t ≥ t 0 , we consider the system x e (t) = (β e − α e )x e − γ e x e + γ c x c , x e (t 0 ) = 0, which has the unique trivial solution solution x c (t) = x e (t) = 0 for t ≥ t 0 , so extinction has occurred at t 0 .It follows that x e x c ≥ 0 on [0, T ], for each 0 < T and so both functions are nonnegative, as long as they exist.
Next, to obtain upper bounds on x c and x e , we add the first two equations and obtain Since x c and x e are nonnegative, each one is bounded and this establishes the first two estimates in (3.6).
We show now that y r is actually positive.We note that y r = γy r (λ − y r ) − δ r x c y r .
Then, since y r ≥ −(γy r + δ r x c )y r , it follows that y r ≥ −(γy r0 e γλt + δ r x co e βct )y r .
This proves both sides of the estimate for y r in (3.6).Finally, we address the last inequality for w.We have Since w 0 > 0, assume that there is 0 < t 0 such that w(t) > 0 for 0 ≤ t < t 0 , and if w(t 0 ) = 0 then η(t 0 ) = 0, and so w (t 0 ) = δ w x c y r ≥ 0. But y > 0 so if x c (t 0 ) = 0, then, as discussed above, both x c and x e vanish identically for t 0 ≤ t, which shows that w(t) = 0 for t 0 ≤ t, as well.Furthermore, it follows from the equation that w ≤ δ w x c y r , and using the estimates for y r and x c , we obtain w ≤ δ w y r0 (x c0 + x e0 )e (γλ+βc)t .
Integration over [0, t] yields w(t) ≤ (w 0 + M ) + M e (γλ+βc)t , where We conclude that w is nonnegative and bounded.This completes the proof.
These estimates allow us to establish the existence of solutions of the system by using Theorem 3.1.Indeed, the a-priori estimates, established in Proposition 3.2, for possible solutions of Model 2.1, allow us to apply Theorem 3.1, hence we obtain the existence of the solutions to the HANDY-SMRN Model.Moreover, a careful checking of the proof of Proposition 3.2 shows that the results hold true even when x 0e = 0 and w 0 = 0, by allowing x 0e → 0 and w 0 → 0. We summarize these results in our main existence theorem.Theorem 3.3.Assume that x c0 and y r0 are positive and x e0 , y n0 and w 0 are nonnegative.Then, Model 2.1 has solutions on every finite time interval [0, T ].Moreover, the solutions satisfy estimates (3.6).
As was mentioned above, Theorem 3.1 does not require the right-hand side F to be Lipschitz continuous, but then, Theorem 3.3 does not guarantee the uniqueness of the solution.That has to be shown separately, and is left open in this work.

Steady states and stability
This section studies the steady states of system (2.1).First, we determine the steady states, and then study their stability.4.1.Steady states.Clearly, the origin z 0 = (0, 0, 0, 0, 0) in R 5  + and z λ = (0, 0, λ, 0, 0) are steady states.As we show below, the origin is unstable while z λ is stable.Then, we turn to find the other possible steady states z = (x c , x e , y r , y n , w), and assume that y r > 0, since without resources the model does not make sense, and xc > 0, since without workers there is no wealth generation.
To find the steady states, we solve the nonlinear algebraic system where C c and C e are given in (2.9), w th is given in (2.7), and α c and α e are given in (2.10).Next, we omit the in w th , since xc > 0, and use We introduce two nonnegative parameters that simplify the writing in what follows: which measure the ratio of wealth to the wealth threshold, and the ratio of Elites to Commoners.Then, C c = min(1, η)sx c , C e = min(1, η)κ 1 sx e , and , We recall that, by assumption, 0 We note that κ 1 ≥ 1 and consider three cases: (i) η < 1/κ 1 , in which there is not enough wealth for the minimal existence of both populations; (ii) 1/κ 1 ≤ η < 1, in which there is enough wealth for the rich but not enough for the workers; and (iii) 1 ≤ η in which there is enough wealth for both populations.
Case (i): η < 1/κ 1 -not enough wealth.Then, C c = ηsx c , C e = ηκ 1 sx e , and Since x c > 0 by assumption, the first equation in (4.1) implies, From the second equation, we get either x e = 0, i.e., ψ = 0 (no Elites), or Adding the two expressions and using simple algebra yields (ψ > 0), It follows that when this inequality does not hold, Case (i) cannot be realized as a nonnegative steady solution.Furthermore, we find that the ratio ψ is determined by the coefficients We note in passing that ψ = 0, i.e., no Elites, can happen only when γ c = 0, so no upward mobility, and β e = β c = α M .We describe this case below.
The equation for y r leads to and since we assume that y r > 0, we find Hence, The equation for w yields and since x c > 0, 0 = δ w y r − ηs − ηκ 1 sψ. Therefore, This implies that Next, x e = ψx c .Finally, w = ηρ(x c + κ 0 x e ), thus Moreover, since x c > 0, we find that w > 0 and which provides the following bound on ψ, Case (ii): 1/κ 1 ≤ η < 1 not enough wealth for the workers, but enough for the rich.Proceeding as above, we find Similarly, since x c > 0, we obtain from the first equation in (2.1), The second equation yields either x e = 0, which is the case with ψ = 0, or 0 = (β e − α m − γ e )ψ + γ c .
Thus, ψ is also fixed in this case, Therefore, 0 ≤ ψ only when 0 < α m − β e + γ e .Otherwise, Case (ii) does not yield any steady states.Then, η is determined as,

.8)
In this case, the right-hand side must satisfy 1/κ 1 ≤ η < 1, otherwise Case (ii) cannot have steady solutions.Proceeding as in Case (i), the equations for y r and w lead to the same expressions but with the new η and ψ.We summarize the findings in Cases (i) and (ii) as follows.
Case (iii): 1 ≤ η there is enough wealth for everyone.We note that in this case the equations for x c and x e are not coupled to those for y and w.Proceeding as above, we obtain Since x c > 0, the first equation in (2.1) leads to From the second equation, we get either x e = 0, which is the case with ψ = 0, or 0 = (β e − α m − γ e )ψ + γ c .
Hence, Thus, Finally, the equation for w yields 0 = δ w y r − s − κ 1 sψ. Therefore, This implies that Next, x e = ψx c .Finally, w = ηρ(x c + κ 0 x e ), thus Now, by assumption x c > 0, hence w > 0 and which provides the following bound on ψ, We summarize Case (iii) as follows.
Remark 4.3.We note that when 1 ≤ η only the wealth w depends on η, while ψ is fixed by the system parameters.Thus, (x c , x e , y r , y n ) are uniquely determined, so they have the same values, independently of the wealth, as long as 1 ≤ η, i.e., there is sufficient wealth to have the populations be above the famine threshold.

4.2.
Stability.We turn to study the stability of the steady states.To that end, we use the Jacobian matrix J of F ( constructed in the Appendix) evaluated at the steady states.It is found that we have five versions of J corresponding to the two 'simple' steady states z = (0, 0, 0, 0, 0) and z = (0, 0, λ, 0, 0), and three different cases of the steady states derived above.To simplify the expressions, we use the notation In the 'simple' steady states, we have α c = α e = α M , hence ∆α = 0; η = 0 since w = 0; and ψ = 1, hence ψ 0 = 1/(1 + κ 0 ) and ψ 1 = 1 + κ 1 , then it is found (Appendix) that when we substitute z = (0, 0, 0, 0, 0) into J (i) , we obtain 12) The eigenvalues are We note that λ 3 is positive, hence the origin is unstable, actually, repelling.
Next, a similar calculation shows that Then, the eigenvalues λ 1,2 are as above and because of the assumption that 0 ≤ ψ only when 0 < α m − β e + γ e and 0 < α m − β c + γ c , they have negative real parts.Also, λ 4 and λ 5 are the same as above while λ 3 = −γλ < 0. Therefore, the steady state is stable and attracting.
These results show that the model allows for the collapse of the society, so that only nature remains, while the populations decay to zero.
We turn to the other steady states.To simplify the presentation, we use the notation In Case (i) we find Where, J 22 = (β e −α e −γ e )−ηκ 0 κ 1 ∆αψ 0 ψ, J 33 = − γηs δw ψ 1 , and J 53 = γδw δr λ − ηs δw ψ 1 .Similarly, in Case (ii), we have Where, Finally, in Case (iii), we obtain .16) Clearly, the study of the states stability can be done only in particular cases when the constants are given.
Thus, we have closed form expressions for the Jacobian matrices in the five cases.In each of the simulations below, we also included the information about which case is simulated and its stability, when it converges to a steady state.
We note that in the baseline simulations of Case (iii), with a modified death rate so that η ≥ 1, we found This means that the steady state is unstable, and further study may show periodic solutions, non-periodic oscillations or chaotic oscillations.

Simulations
This section presents and briefly discusses the results of our model numerical experiments.The main interest is to exhibit the four different types of behavior of the solutions: monotone or oscillatory approach to the steady state; periodic solutions; and system collapse.It is very likely that the model can predicts chaotic behavior, however, we did not observe chaos with the system parameters we have experimented with.The description of the code for the simulations can be found in [1].
We start with a short study of the effects of increasing δ n , the the renewables depletion coefficient, and the effects of doubling κ 1 , the pay gap measure.
Table 1 lists the model dependent variables, their initial values, and the values of the parameters.These were chosen in the baseline simulation, which depicts a typical system behavior.To obtain the four different types of behavior, we changed some of the parameters and initial conditions, and the changes are described with each type of behavior.The results are depicted in the Figs.1-6.The numbers in the table include multiples of the renewable resources depletion coefficient, to highlight the over consumption of resources.
In all the simulations, the blue curve depicts the workers x c ; magenta the rich x e ; green and black curves are the renewables y r and nonrenewables y n , respectively; and the society's wealth w is represented by the cyan curve.
We note that here, unlike [5], we use 'time units' instead of 'years,' on the horizontal scale, since our units are quite arbitrary, so the time is not necessarily identifiable as a day, month or year.The main purpose is to depict the different types of solutions.1, and all fit Case (i) of the analysis, in which η ≤ 1/κ 1 = 1, since we use κ 1 = 1. Figure 1 (TL-top left) describes a 'soft landing' [7] in which the system approaches the steady state of (45000, 4000, 60, 0, 5.5) monotonically after a single initial oscillation.From time 500 onward, the system converges monotonically to this steady state.The nonrenewables tend to zero, as expected, while the renewable resources tend to 60, which is lower than λ = 100.The two populations and the wealth converge to their limits smoothly.The depletion factor δ r = 1.3 δ * has a different value than in Table 1.Also, η = 0.67 at the steady state.The two populations grow to high values causing complete depletion of the nonrenewables and almost all of the wealth, thus, leaving the society's populations with only the renewables.The wealth decrease to 5.5 from its full value of 20 at the peek of the population rise.The vanishing of the nonrenewables and a considerable decrease in wealth drag the populations to half of their prosperity peeks.
Figure 1 (TR-top right) depicts convergence to the steady state after a long period of large but damped oscillations, and then after time 575, the oscillations almost vanish and the steady state is rapidly approached.We note that, as expected, the nonrenewables do not oscillate, since these are just consumed.The value of δ r = 4.6 δ * was different than that in Table 1, and we increased the time span for a better depiction of the approach to the steady state.The steady state the system settles to around (27000, 25000, 20, 0, 0).Moreover, η = 0.67 at the steady state, which is Case (i).
Figure 1 (BL bottom left) provides numerical evidence that the system may have periodic solutions in (x c , x e , y r , w), once the initial conditions have been chosen appropriately.However, because of the structure of the system, the nonrenewables y n cannot oscillate, and must decay.Here, δ r = 7.2δ * , and a shift of the initial conditions was made to (750, 750, 96, 0, 0) to generate the periodic solution.Moreover; the following changes were also used, (γ e c , γ c e , δ r , δ n , δ w ) = (0.085, 0.075, 7.2 δ * , 6.15 • 10 −6 , 1.18δ r ).The mathematical proof of the existence of periodic solutions is unavailable, yet, however, these simulations provide a strong indication that they are likely to exist.
Finally, Figure 1 (BR) depicts the system collapse at about time 600.The populations increase considerably from the initial conditions and then the resources, the wealth and the populations collapse.
We note that the populations' values at nature's collapse are of interest.For instance; in Figure 1 (BL) the renewables and wealth collapse around time 210.The populations collapse after another hundred time units, when the commoners population decreases to 850; however, the system starts to recover after another 100 time units, and the renewables increase about time 310 with the commoners working to revive the wealth.In the meantime; when δ r is large, such as in Figure 1 (BR), the excessive consumption of resources doesn't allow the system to recover; keeping η close to zero which is an unstable case for the system.5.2.Dependence on δ n , δ r .We study numerically the system's dependence on the parameters δ n and then δ r , varying the rates of depletion and consumption of the nonrenewables and renewables, respectively.This is the new part of the modified model.When the depletion factor of the nonrenewables is decreased, δ n = 6.15 • 10 −8 , i.e. 100-fold less than δ n = 6.15 • 10 −6 , which means a more moderate use of nonrenewables, we find, Figure 2(L), that the system a steady state after 2000 time units, and the populations' numbers are considerably larger than in the case of increased depletion, δ n = 6.15 • 10 −5 , Figure 2 (R).In both cases the system eventually reaches its steady state, but, in both cases the populations almost collapse following the collapse of the renewables.We mention that a decrease in δ n causes the wealth to increase; that is also the reason behind the expansion of the two populations.That is not the case with an increased δ n , shown in Figure 2 (R), since the nonrenewables are drained because of the overuse in a shorter time period, so the system depends only on the renewables as a source of wealth.In both figures the change was done only to δ n ; while keeping the constants at the values (γ e c , γ c e , δ r , κ 1 , δ w ) = (0.085, 0.075, 2.19 δ * , 1, 1.18δ r ).To obtain some insight into the effects of increasing the use of renewable resources, we performed an experiment by changing the renewables rate constant δ r , and the results are provided in Figures 3-5.Here, we also use κ 1 = 2.It is found that within the whole interval 9.96δ * ≤ δ r ≤ 23.2δ * there are rapid oscillations.However, Figure 3 (L), where δ r = 9.9δ * , seems to be of interest, since there are many fast slightly damped oscillations until about time 7100, and then, in a very short period of about 100 time units, the system approaches the steady state.The solution about time 7000, and the fast approach to the steady state, are depicted in the zoom 3 (R).On the other hand, in Figure 4 we use δ r = 23.2δ* and it seems that the solutions are periodic.We use the same parameters in both Figures 3-5, except for δ r .For the sake of completeness, in a shorter time period of 5000 time, Figure 5 depicts the death and consumption rates as functions of time for the scenario with δ r = 23.2δ* .Moreover, it is found that during the oscillations, the system approaches partial collapse and then it recovers.
We note an interesting behavior in Figure 5 (L), since having piecewise constant values in a dynamical system, for the consumption rates, is somewhat unusual, and here it seems to be related to the non-differentiability of the rate functions.Also, the death rates have rather interesting spikes.5.3.Dependence on κ 1 .Finally, we study the system's dependence on the consumption ratio κ 1 ; noting that in the previous simulations, we used κ 1 = 1.Here, we used four cases κ 1 = 0.5, 1, 5, 10.The simulations results are depicted in Figure 6.We use γ e = 0.075 and γ c = 0.085.The case κ 1 = 0.5 is both unrealistic in terms of the model, since the rich would not consume less than the commoners, and mathematically it produces what seems to be a periodic solution.The other three cases lead to very similar approach to the steady states and comparable steady states, (z = x c , x e , y r , y n , w).The convergence is either monotone or with a few oscillations before the monotone convergence.Indeed, the steady states were found to be (approximately), κ 1 = 1, z = (29000, 33000, 23, 0, 4); κ 1 = 5, z = (27500, 41600, 27, 0, 1); κ 1 = 10, z = (22600, 35000, 40, 1).However, we note that the ways they were reached seem to be quite different, with more or less societal upheaval in the form of large oscillations.This indicates that just studying the steady states of a society's model may not provide a good description and insight into the society's dynamics.

Conclusions and further research
This work presents the HANDY-SMRN model, which is a new version of the HANDY model [5], that is more realistic.It is a continuation of the studies of other versions of the HANDY model in [1].This HANDY-SMRN model extends the HANDY model by taking into account mobility between the Commoners (workers) and Elites (rich), and splits the natural resources into renewables and nonrenewables, where the latter once consumed cannot be renewed.The model is a nonlinear coupled system of five ODEs for the Commoners and Elite populations, societal wealth, and renewable and nonrenewable resources.It establishes the existence, positivity and boundedness of the solutions of the model.We note that similar properties of the HANDY model have been established in [1], and the analysis here is similar, but tailored to the HANDY-SMRN model.Then, it studies the system's steady states and their stability.It is found that the steady states depend essentially on the ratio η = w/w th of the wealth w to the wealth's threshold w th .This type of models, which describe in a very coarse way complex societies, although relatively 'simple,' seem to provide insight into what may be common in the collapse of such societies.Clearly, there is considerable interest to find out if the model can fit some of the historical collapses of societies which essentially consisted of commoners and nobility.
Mathematically, there is interest in establishing that the model has oscillatory solutions (except for the nonrenewables), since our computer experiments indicate that this is the case.Moreover, it is very likely that a fifth type of behavior exists, namely chaotic behavior, and it is of interest to obtain numerically chaotic solutions, and possibly to prove the existence of chaos rigorously.However, these questions remain unresolved, yet, and will need future investigations.
The area is wide open and there is plenty of room for various types of mathematics that may shed light on social dynamics.

Figure 4 .
Figure 4.A scenario with δ r = 23.2δ* ; (L) the system settles quickly into periodic oscillations with period of about 1000 time units; (R) zoom detail.

C
c = ηsx c , C e = ηκ 1 sx e , α c = α M − η(α M − α m ), α e = α m .It follows, since α e = α m , that the only changes from Case (i) are those related to α e .Thus, we find that the only differences in (4.14) are in the following components:J 21 = γ c , J 22 = β e − α m − γ e , J 25 = 0,and this yields (4.15).Finally, in Case (iii), we have C c = sx c , C e = κ 1 sx e , α c = α e = α m .Then, the changes in J(z), as compared to Case (i) are those related to α c and α e .Therefore,

Table 1 .
Parameters and data used in the simulations of the HANDY-SMRN model.Qualitative behavior.The four model qualitatively different types of solutions are depicted in Figure