Lotka-Volterra Equations Revisited

In this section, we will revisit the Lotka-Volterra Systems discussed in the first chapter. However, this time we will create system models from individual components. After recreating the behavior shown in the first chapter, we’ll expand the set of effects we consider and reconfigure these component models into other system models that demonstrate different dynamics.

Classic Lotka-Volterra

We’ll start by looking at the classic Lotka-Volterra system. In order to create such a system using component models, we will require models to represent the population of both rabbits and foxes as well as models for reproduction, starvation and predation.

Connectors

However, as we learned during our discussion of Connectors, before we can start building component models we first need to formally define the information that will be exchanged by interacting components by defining connectors. The connector we will use in this section is the Species connector and it is defined as follows:

within ModelicaByExample.Components.LotkaVolterra.Interfaces;
connector Species "Used to represent the population of a specific species"
  Real population "Animal population";
  flow Real rate "Flows that affect animal population";
end Species;

This connector definition is interesting because these definitions do not come from engineering. Instead, they really arise from ecology. In this case, our across variable is population which represents the actual number of animals of a particular species. Our through variable, indicated by the presence of the flow qualifier, is rate which represents the rate at which new animals “enter” the component that this connector is attached to.

Regional Population

To track the population of a given species in one region, we’ll use the RegionalPopulation model. The model has several noteworthy aspects so we’ll present the model piece by piece starting with:

within ModelicaByExample.Components.LotkaVolterra.Components;
model RegionalPopulation "Population of animals in a specific region"
  encapsulated type InitializationOptions = enumeration(
      Free "No initial conditions",
      FixedPopulation "Specify initial population",
      SteadyState "Population initially in steady state");

The first two lines are as expected. But after that we see that this model defines a type called InitializationOptions. The type definition is qualified with the encapsulated keyword. This is necessary because the type is being defined within a model and not a package. Modelica has a rule that if we wish to refer to this type from outside the model definition, the type definition must be prefixed by encapsulated. We can see from the definition of this enumeration that it defines three distinct values: Free, FixedPopulation and SteadyState. We’ll see how these values will be used shortly.

The first declarations in our RegionalPopulation model is:

  parameter InitializationOptions init=InitializationOptions.Free
    annotation(choicesAllMatching=true);

Note that the first parameter, init, utilizes the InitializationOptions enumeration both to specify its type (the enumeration itself) and its initial value, Free. Also note the presence of the choicesAllMatching annotation. We’ll talk more about this annotation later in this chapter, when we are reviewing the concepts introduced here, and in subsequent chapters.

The next declaration is:

  parameter Real initial_population
    annotation(Dialog(group="Initialization",
    enable=init==InitializationOptions.FixedPopulation));

The initial_population parameter is used to represent the initial value for the population in this region at the start of the simulation. However, as we will see shortly in the equation section, this value is only used if the value of init is set to FixedPopulation. For this reason, the enable annotation on this parameter is set to init==InitializationOptions.FixedPopulation. This annotation is used to inform Modelica tools of this relationship. This information can then be taken into account when building graphical user interfaces (e.g., a parameter dialog) associated with this model.

It is also worth noting the presence of the Dialog annotation in the definition of initial_population. This annotation allows the model developer to organize parameters into categories, in this case “Initialization”. Tools generally use such information to help structure parameter dialogs.

The last public declaration in the model is for the connector instance that allows interactions with other components:

  Interfaces.Species species
    annotation (Placement(transformation(extent={{-10,90},{10,110}}),
        iconTransformation(extent={{-10,90},{10,110}})));

Here we again see the Component Placement annotation which we will again defer talking about for the moment. This leaves the last declaration, which happens to be protected:

protected
  Real population(start=10) = species.population "Population in this region";

This variable represents the number of animals in this region. It is given a non-zero start value to avoid the trivial solutions we saw in our earlier discussion of Steady State Initialization of Lotka-Volterra systems. We can also see, from this declaration, that this declaration equates the local variable, population, with the value of the across variable on the species connector, species.population. In effect, the population variable is an alias for the expression species.population.

Now that we have the declarations out of the way, let’s look at the equations associated with the RegionalPopulation model:

initial equation
  if init==InitializationOptions.FixedPopulation then
    population = initial_population;
  elseif init==InitializationOptions.SteadyState then
    der(population) = 0;
  else
  end if;
equation
  der(population) = species.rate;
  assert(population>=0, "Population must be greater than zero");
end RegionalPopulation;

The initial equation demonstrates the significance that the value of the init parameter has on the behavior of this component. In the case that the init value is equal to the FixedPopulation value in the InitializationOptions enumeration, an equation is introduced specifying that the value of the population variable at the start of the simulation is equal to the initial_population parameter. If, on the other hand, the value of init is equal to the SteadyState value of the enumeration, then an equation is introduced specifying that the rate of population change at the start of the simulation must be zero. Obviously, if init is equal to Free (the last remaining possibility), no initial equation is introduced.

Within the equation section, we see that the rate at which the population changes is equal to the value of the flow variable on the species connector, species.rate. Again, recall the sign convention that a positive value for a flow variable means a flow into the component and the fact that this equation is consistent with that sign convention (i.e., a positive value for species.rate will have the effect of increasing the population within the region).

The last thing worth noting about the model is the presence of the assert call in the equation section. This is used to define a model boundary (i.e., a point beyond which the equations in the model are not valid). It is used to enforce the constraint that the population within a region cannot be less than zero. If a solution is ever found where the constraint is violated, the resulting error message will be “Population must be greater than zero”.

This model also has an Icon annotation associated with the model definition. As usual, the Icon annotation is not included in the source listing. But when this component model is rendered within a system model, its icon will look like this:

Reproduction

The first real effect we will examine is reproduction. As we know from our previous discussion, the growth in a given population due to reproduction is proportional to the number of animals of that species in a given region. As a result, we can describe reproduction very succinctly as:

within ModelicaByExample.Components.LotkaVolterra.Components;
model Reproduction "Model of reproduction"
  extends Interfaces.SinkOrSource;
  parameter Real alpha "Birth rate proportionality constant";
equation
  growth = alpha*species.population "Growth is proporational to population";
end Reproduction;

where alpha is the proportionality constant. However, the simplicity and clarity of this model is due mostly to the inheritance of the SinkOrSource model in much the same way that our “DRY” Electrical Components benefited from inheriting the TwoPin model.

The SinkOrSource model is a starting point for any model that either creates or destroys animals in a population. It is defined as follows:

within ModelicaByExample.Components.LotkaVolterra.Interfaces;
partial model SinkOrSource "Used to describe single species effects"

  Species species
    annotation (Placement(transformation(extent={{-10,90},{10,110}})));
protected
  Real growth "Growth in the population (if positive)";
  Real decline "Decline in the population (if positive)";
equation
  decline = -growth;
  species.rate = decline;
end SinkOrSource;

To understand these equations it is first necessary to understand that any model that extends from SinkOrSource will generally be connected to a RegionalPopulation instance (but will not, itself, be a RegionalPopulation model). This means that if the flow variable species.rate in such an instance is positive, it will have the effect of pulling animals out of the RegionalPopulation model. Looking at the SinkOrSource model in this way, we can see that the variable decline is simply an alias for species.rate. In other words, when decline has a positive value, species.rate will have a positive value and, therefore, any RegionalPopulation that this SinkOrSource instance is connected to will suffer a drain on its population. Conversely, the growth variable is positive when species.rate is negative. In that case, the connected RegionalPopulation model will see an increase in species population.

By defining the SinkOrSource model and inheriting from it, much of this complexity is hidden. As a result, models like Reproduction can have equations written in a way that make their behavior more intuitive, e.g., growth = alpha*species.population.

Although not shown, the Icon for the Reproduction model is rendered as:

Starvation

Just like the Reproduction model just described, the Starvation model also inherits from the SinkOrSource model. However, its behavior with respect to the decline variable, is described as follows:

within ModelicaByExample.Components.LotkaVolterra.Components;
model Starvation "Model of starvation"
  extends Interfaces.SinkOrSource;
  parameter Real gamma "Starvation coefficient";
equation
  decline = gamma*species.population
    "Decline is proporational to population (competition)";
end Starvation;

Predation

The last effect we need to consider before building a system model to represent the classic Lotka-Volterra behavior is a model for predation.

Recall our previous discussion of the SinkOrSource model and the potential confusion associated with sign conventions. The SinkOrSource model was designed to work with effects that only interacted with a single RegionalPopulation (since it had only one Species connector). In order to address the same potential sign convention confusion for effects that involve interactions between two different regional populations, the following partial model, Interaction was defined:

within ModelicaByExample.Components.LotkaVolterra.Interfaces;
partial model Interaction "Used to describe interactions between two species"
  Species a "Species A"
    annotation (Placement(transformation(extent={{-110,-10},{-90,10}})));
  Species b "Species B"
    annotation (Placement(transformation(extent={{90,-10},{110,10}})));
protected
  Real a_growth "Growth in population of species A (if positive)";
  Real a_decline "Decline in population of species A (if positive)";
  Real b_growth "Growth in population of species B (if positive)";
  Real b_decline "Decline in population of species B (if positive)";
equation
  a_decline = -a_growth;
  a.rate = a_decline;
  b_decline = -b_growth;
  b.rate = b_decline;
end Interaction;

Again, we have the concepts of growth and decline variables. However this time we have two version of each. One is associated with the a connector and the other is associated with the b connector.

Using these definitions, we can define Predation very simply as:

within ModelicaByExample.Components.LotkaVolterra.Components;
model Predation "Model of predation"
  extends Interfaces.Interaction;
  parameter Real beta "Prey (Species A) consumed";
  parameter Real delta "Predators (Species B) fed";
equation
  b_growth = delta*a.population*b.population;
  a_decline = beta*a.population*b.population;
end Predation;

This model captures the effect that the growth in the “B” (predator) population is proportional to the product of the predator and prey populations. Similarly, the decline in the “A” (prey) population is also proportional the product of the predator and prey populations (although with a different proportionality constant).

Although not shown, the Icon for the Predation model is rendered as:

Note that the Predation model is asymmetric. The b connector should be connected to the predator population and the a connector should be connected to the prey population. This is reinforced by the image and asymmetry of the icon itself.

Classic System Model

With all of these components in hand, we can very easily construct a component-oriented version of the classic Lotka-Volterra behavior by dragging and dropping the components into the following system configuration:

Component-oriented version of the classic Lotka-Volterra model

Here we see that the Starvation model is attached to the foxes population while the Reproduction model is attached to the rabbits population. The Predation model is connected to both populations with the a (prey) connector attached to the rabbits and the b (predator) connector attached to the foxes.

As we can see from the following plot, the behavior of this system is identical to the one presented in our earlier discussion of Lotka-Volterra Systems:

../../../_images/CLV.png

Introducing a Third Species

As we will see over and over again, there is an initial investment in building component models over simply typing the equations that they encapsulate. But there is also a significant “payoff” to this process because of the schematic based system composition that is then possible as a result. In the context of the Lotka-Volterra example, this is exemplified by the addition of a third species, wolves, to the classic Lotka-Volterra system.

Adding Wolves

The creation of a model with a third species does not require any additional component models to be defined. Instead, we can reuse not only our existing models for Predation, Starvation and RegionalPopulation, but we can also reuse the ClassicLotkaVolterra model itself:

within ModelicaByExample.Components.LotkaVolterra.Examples;
model ThirdSpecies "Adding a third species to Lotka-Volterra"
  import ModelicaByExample.Components.LotkaVolterra.Components.RegionalPopulation.InitializationOptions.FixedPopulation;
  extends ClassicLotkaVolterra(rabbits(initial_population=25), foxes(initial_population=2));
  Components.RegionalPopulation wolves(init=FixedPopulation, initial_population=4)
    annotation (Placement(transformation(extent={{30,40},{50,60}})));
  Components.Starvation wolf_starvation(gamma=0.4)
    annotation (Placement(transformation(extent={{70,40},{90,60}})));
  Components.Predation wolf_predation(beta=0.04, delta=0.08) "Wolves eating Foxes"
    annotation (Placement(transformation(extent={{-10,-10},{10,10}}, rotation=90, origin={60,20})));
  Components.Predation wolf_rabbit_predation(beta=0.02, delta=0.01) "Wolves eating rabbits"
    annotation (Placement(transformation(extent={{-10,30},{10,50}})));
equation
  connect(wolf_predation.b, wolves.species) annotation (Line(
      points={{60,30},{60,80},{40,80},{40,60},{40,60}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolf_rabbit_predation.a, rabbits.species) annotation (Line(
      points={{-10,40},{-40,40},{-40,-20}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolf_predation.a, foxes.species) annotation (Line(
      points={{60,10},{60,0},{40,0},{40,-20}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolf_starvation.species, wolves.species) annotation (Line(
      points={{80,60},{80,80},{40,80},{40,60}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolves.species, wolf_rabbit_predation.b) annotation (Line(
      points={{40,60},{40,80},{20,80},{20,40},{10,40}},
      color={0,127,0},
      smooth=Smooth.None));
  annotation (experiment(StopTime=100, Tolerance=1e-006));
end ThirdSpecies;

Such a model would not typically be created by typing in the source code you see above. Instead, within a graphical development environment it would take less than a minute to assemble such a system by augmenting the existing ClassicLotkaVolterra model. When visualized, the schematic for the resulting system is rendered as:

The Classic Lotka-Volterra model augmented with an additional predatory species

Resulting Dynamics

By creating such a model, we can quickly explore the differences in system dynamics between the classic model and this three species model. The following plot shows how these three species interact:

../../../_images/ThirdS.png

By using the init parameter in the various RegionalPopulation instances, we can also quickly create a model to solve for the equilibrium population levels for all three species:

within ModelicaByExample.Components.LotkaVolterra.Examples;
model ThreeSpecies_Quiescent "Three species in a quiescent state"
  import ModelicaByExample.Components.LotkaVolterra.Components.RegionalPopulation.InitializationOptions.SteadyState;
  extends ThirdSpecies(
    rabbits(init=SteadyState, population(start=30)),
    foxes(init=SteadyState, population(start=2)),
    wolves(init=SteadyState, population(start=4)));
  annotation (experiment(StopTime=100, Tolerance=1e-006));
end ThreeSpecies_Quiescent;

All that is required in this model is to extend from the ThirdSpecies model and modify the init parameter for each of the underlying species populations. Simulating this model gives us the equilibrium population level for each species:

../../../_images/ThreeS.png

From an ecological standpoint, we can already make an interesting observation about this system. If we start it from a non-equilibrium condition the system quickly goes unstable. In other words, the introduction of wolves into the otherwise stable eco-system involving only rabbits and foxes has a significant impact on the population dynamics.