# 4.6: Building Your Own Model Equations with Multiple Variables

We can take one more step to increase the complexity of the model building, by including more than one variable. Following the theme of population growth, let’s consider ecological interactions between two species. A typical scenario would be the *predator-prey interaction*. Let’s stick to the population-level description of the system so each species can be described by one variable (say, \(x\) and \(y\)). The ﬁrst thing you should consider is each variable’s inherent dynamics, i.e., what would happen if there were no inﬂuences coming from other variables. If there is always plenty of food available for the prey, we can assume the following:

- Prey grows if there are no predators.

- Predators die if there are no prey.

These assumptions can be diagrammatically illustrated in Fig. 4.6. This type of diagram is called a *causal loop diagram* in* System Dynamics* [26]. Each circle, or node, in this diagram represents a state variable in the system. The self-loop arrow attached to each node represents the effect of the variables on itself (e.g., the more prey there are, the faster their growth will be, etc.). The plus/minus signs next to the arrows show whether the effect is positive or negative. We can now consider the interactions between the two variables, i.e., how one variable inﬂuences the other and vice versa. Naturally, there should be the following effects:

- The prey’s death rate increases as the predator population increases.
- The predators’ growth rate increases as the prey population increases.

These interactions can be illustrated as arrows between nodes in the causal loop diagram (Fig. 4.7).

Now is the time to translate the structure of the system illustrated in the diagram into mathematical equations. Each arrow in the diagram tells us whether the effect is positive or negative, but they don’t give any exact mathematical form, so we will need to create a mathematical representation for each (possibly using the aforementioned tips).

The inherent dynamics of the two variables are quite straightforward to model. Since we already know how to model growth and decay, we can just borrow those existing models as building components, like this:

\[x_{t} =x_{t-1}+r_{x}x_{t-1}(1-\frac{x_{t-1}}{K})\label{(4.28)}\]

\[y_{t}=y_{t-1} -d_{y}y_{t-1}\label{(4.29)}\]

Here, I used the logistic growth model for the prey \((x)\) while using the exponential decay model for the predators \((y)\). \(r_{x}\) is the growth rate of the prey, and \(d_{y}\) is the death rate of the predators \((0 < dy < 1)\).

To implement additional assumptions about the predator-prey interactions, we need to ﬁgure out which part of the equations should be modiﬁed. In this example it is obvious, because we already know that the interactions should change the death rate of the prey and the growth rate of the predators. These terms are not yet present in the equations above, so we can simply add a new unknown term to each equation:\[x_{t} =x_{t-1} +rx_{t-1}(1-\frac{x_{t-1}}{K})-d_{x}(y_{t-1})x_{t-1}\label{(4.30)}\]

\[y_{t}=y_{t-1} - dy_{t-1} +r_{y}(x_{t-1})y_{t-1}\label{(4.31)}\]

Now the problems are much better deﬁned. We just need to come up with a mathematical form for \(d_{x}\) and \(r_{y}\).

The death rate of the prey should be 0 if there are no predators, while it should approach 1 (= 100% mortality rate!) if there are too many predators. There are a number of mathematical formulas that behave this way. A simple example would be the following hyperbolic function

\[d_{x}(y)=1-\frac{1}{by+1},\label{(4.32)}\]

where \(b\) determines how quickly \(d_{x}\) increases as \(y\) increases. The growth rate of the predators should be 0 if there are no prey, while it can go up indeﬁnitely as the prey population increases. Therefore, the simplest possible mathematical form could be

\[r_{y}(x) =cx,\label{(4.33)}\]

where \(c\) determines how quickly \(r_{y}\) increases as \(x\) increases. Let’s put these functions back into the equations. We obtain the following:

\[x_{t} = x_{t-1} +rx_{t-1}(1-\frac{x_{t-1}}{K})-(1-\frac{1}{by_{t-1} +1})x_{t-1}\label{(4.34)}\]

\[y_{t} =y_{t-1}-dy_{t-1} +cx_{t-1}y_{t-1}\label{(4.35)}\]

Exercise 4.11

Test the equations above by assuming extreme values for \(x\) and \(y\), and make sure the model behaves as we intended.

Exercise 4.12

Implement a simulation code for the equations above, and observe the model behavior for various parameter settings.

Figure 4.8 shows a sample simulation result with \(r = b = d = c = 1\), \(K = 5\) and \(x_{0} = y_{0} = 1\). The system shows an oscillatory behavior, but as its phase space visualization (Fig. 4.8 right) indicates, it is not a harmonic oscillation seen in linear systems, but it is a nonlinear oscillation with distorted orbits.

The model we have created above is actually a variation of the *Lotka-Volterra model*, which describes various forms of predator-prey interactions. The Lotka-Volterra model is probably one of the most famous mathematical models of nonlinear dynamical systems that involves multiple variables.

Exercise 4.13

Try several different parameter settings for \(r,b,d\), and \(c\), and see how the system’s behavior changes. In some cases you may ﬁnd an unrealistic, invalid behavior (e.g., indeﬁnite growth of predators). If so, revise the model to ﬁx the problem.

In this chapter , we reviewed some basics of mathematical modeling in difference equations. As I keep saying, the best way to learn modeling is through practice. Here are some more modeling exercises. I hope you enjoy them!

Exercise 4.14

Develop a discrete-time mathematical model of two species competing for the same resource, and simulate its behavior.

Exercise 4.15

Consider the dynamics of public opinions about political ideologies. For simplicity, let’s assume that there are only three options: conservative, liberal, and neutral. Conservative and liberal are equally attractive (or annoying, maybe) to people, with no fundamental asymmetry between them. The popularities of conservative and liberal ideologies can be represented by two variables, \(p_{c}\) and \(p_{t}\), respectively \((0 ≤ pc ≤ 1; 0 ≤ pl ≤ 1; 0 ≤ pc + pl ≤ 1)\). This implies that \(1−p_{c} −p_{l} = p_{n}\) represents the popularity of neutral.

Assume that at each election poll, people will change their ideological states among the three options according to the irrelative popularities in the previous poll. For example, the rate of switching from option \(X\) to option \(Y\) can be considered proportional to \((p_{Y} − p_{X})\) if \(p_{Y} > p_{X}\), or 0 otherwise. You should consider six different cases of such switching behaviors (conservative to liberal, conservative to neutral, liberal to conservative, liberal to neutral, neutral to conservative, and neutral to liberal) and represent them in dynamical equations. Complete a discrete-time mathematical model that describes this system, and simulate its behavior. See what the possible ﬁnal outcomes are after a sufﬁciently long time period.

Complete a discrete-time mathematical model that describes this system, and simulate its behavior. See what the possible ﬁnal outcomes are after a sufﬁciently long time period.

Exercise \(4.16\)

Revise the model of public opinion dynamics developed in the previous exercise so that the political parties of the two ideologies (conservative and liberal) run a political campaign to promote voters’ switching to their ideologies from their competitions, at a rate \(inversely\) proportional to their current popularity (i.e., the less popular they are, the more intense their campaign will be). Simulate the behavior of this revised model and see how such political campaigning changes the dynamics of the system.