Assignment 1: Simulating Ecological Homeostasis¶
Background: The Gaia Hypothesis & Daisy World¶
In the 1970s, James Lovelock and Lynn Margulis developed the Gaia Hypothesis. This is a highly controversial theory which proposes that the global ‘biosphere’ is self-regulating— that ecosystems work to maintain the conditions necessary life to persist on Earth. As a case in point to support his theory, Lovelock observed that the amount of heat arriving from our Sun has increased significantly since the origin of life, but that the tempera- ture of the planet has remained within the limits of hability throughout. How does this homeostasis (loosely translatable as ‘keeping things the same’) occur?
Lovelock argues that the homeostasis is a natural biproduct of complex ecosystems, but this is a very controversial claim, because it seems to contradict classic ideas of evolution being driven by selfish competition and ‘survival of the fittest.’ Lovelock’s theory suggests that instead of competing and ‘being selfish,’ diverse species come to some sort of agreement, whereby they’ll cooperate to make conditions beneficial for all life (and not just themselves). If you are interested to know more, the wikipedia page on the ‘Gaia Hypothesis’ provides some more information including links to further reading.
Debate concerning the Gaia Hypothesis continues today. A computational model known as ‘Daisy World’ has been developed to show how regulation of environmental conditions (such as the temperature of the planet) can happen. The model considers a planet populated by two different types of daisies, black ones and white ones. The black daisies increase the temperature (by a small amount) by absorbing sunlight, and the white ones reflect light back into space, decreasing the local temperature. In this assignment, we will develop a simplified version of Daisy World that only includes the white daisies.
As is common in developing computational models, we will make many simplifying assumptions. We will assume:
1. the daisies grow better at certain temperatures than others;
2. the temperature of the planet is the same everywhere;
3. daisies are the only relevant species that lives on the planet;
4. the carrying capacity (i. e. the number of daisies that are alive at equilibrium) is well approximated as a function of temperature;
5. when there are no daisies around, it is still possible for a population to grow (we can imagine there are always seeds in the ground).
Model¶
Our model will be a continuous-time dynamical system that is described by three differ- ential equations. Each equation describes the dynamics of one of the following variables.
• $n$ — the normalized population density of daises. 1 means the planet is covered in the maximum possible density of daises and 0 means there are no daisies at all, etc.
• $y$ – the temperature of the planet (this corresponds in our simplified model to the temperature everywhere)
• $p$ – an external ‘perturbing’ force to the temperature. This can be conceived of as an increase in the radiation from the star that the planet is orbiting, or any other (combination) of factors that are modifying the planet’s temperature other than the presence of the daisies.
In this assignment, as you follow the instructions, you will build a computational model of this system. By the end of it hopefully you will have gained some insight into one way that biological-feedback can (at least in theory) stabilise environments as claimed by Lovelock and Margulis.
Step 1: Write a function that maps the temperature to the carrying capacity of the daisies.¶
In our equations, the carrying capacity is written $k$. It describes the equilibrium population density, i. e. what the population density (n) would move toward if everything else were held constant. Assumption 4 tells us that we can $k$ is well approximated as some function of temperature ($y$):
$$ k = f(y) $$
A typical response of a carrying capacity to an environmental feature such as temperature is a Gaussian function:
$$ k = f(y) = a \cdot \exp \left ( -\frac{(y-b)^2}{2c^2} \right) $$
In this equation, there are three parameters ($a$, $b$ and $c$). $a$ specifies the height of the curve’s peak; $b$ specifies the position of the maximum; and $c$ specifies the standard deviation, i. e. the ‘width’ of the curve.
What values for these parameters should we use? $n$ is normalized, so the peak of the curve should be 1. Let us assume that the daisies grow best at a temperature of $25^o$ C. and that the standard deviation is 3.
Q1a. Write a Python function that takes y as an argument and returns n.
In [1]:
## Q1a answer goes here…
Q1b. To confirm that your function is correct, plot it for suitable values of y.
FOR ALL PLOTS LABEL AXES AND PROVIDE TITLES OR POINTS WILL BE LOST!
In [3]:
## Q1b answer goes here…
Step 2: Describing how daisy density changes over time.¶
In our model, the population density ($n$) changes as a function of the current population density ($n$) and the current carrying capacity ($k$):
$$ \frac{dn}{dt} = g(n,k) $$
We shall model $n$ as always approaching the carrying capacity, with a rate of this approach that is linearly proportional to the difference between the current population density and the carrying capacity. The scale of the response (i. e., the slope of the line), $s = 0.5$.
Q2a. Write down (in math notation, not code) the differential equation that captures what has just been described above. If it is not clear to you how to do this, think about what a graph of this function should look like. Use pencil and paper to think this through.
Double click on the word ANSWER below, to edit the markdown. Your answer here should sit between the two sets of \$\$ \$\$ to be displayed as latex math. When you’ve written an equation, to see it rendered as math, run the cell (I hit control-enter to do so).
The basics of latex math in markdown can be found here: https://www.latex-tutorial.com/tutorials/amsmath/
$$ \frac{dn}{dt} = g(n,k) = ANSWER $$
Q2b. Translate the above mathematical equation into a Python function that takes two arguments (n and k) and returns the rate at which n changes. In other words, fill in the following blank.
In [ ]:
## Q2b answer goes here…
def dndt(n,k):
return …
Q2c. To confirm that your function is working how it should be, plot $g (n,k)$ for an appropriate range of values of n, with k = 0.5.
In [ ]:
## Q2c answer goes here…
Q2d. Super-impose on the same plot g ( n, k ) when k is fixed to one or two other values (you pick the values).
In [ ]:
## Q2d answer goes here…
Q2e. In a sentence or two, describe why the above plots confirm that this function is working correctly. DON’T FORGET TO LABEL YOUR PLOTS.
Step 3: Describing how temperature changes over time.¶
The temperature changes at a rate proportional to the sum of (1) the influence of every-thing other than the daisies, which is captured by the variable $p$, and (2) the influence of the daisies. The influence of the daisies will be linearly proportional to their density, giving us the following differential equation,
$$ \frac{dy}{dt} = p + rn $$
where $r$ is a constant that describes the amount of influence that the daisies have. We’ll assign this parameter a value of $r = -10$, which can be interpreted as saying that when the daisies are at their maximum possible density, they bring the global temperature down by 10 degrees.
Q3. Translate the above Equation into a Python function as you did with the differential equation for $n$.
In [7]:
## Q3 answer goes here…
Step 4: Non-daisy influence¶
Let us imagine that the radiation of the sun is increasing at an increasing rate. This can be described using the following differential equation, where $q = 0.2$ is a constant representing the rate of that increase and $t$ represents the amount of time that has passed since the start of the simulation.
$$ \frac{dp}{dt} = h(t) = qt $$
Q4. Translate this differential Equation into a Python function as you have done with the differential equations for $n$ and $y$.
In [9]:
## Q4 answer goes here…
At this point you should have four Python functions:
def k(…):
“””returns carrying capacity”””
…
def dndt(…):
“””returns rate of change of population density “””
…
def dydt(…):
“””returns rate of temperature change”””
…
def dpdt(…):
“””returns rate of change of non-daisy infl. upon temperature”””
….
Q5. Write a trajectory function that uses these functions to simulate a trajectory from the initial condition that is passed in as an argument.
Your function should use the Euler integration we have learned about in class, and must NOT use an existing library (e.g. scipy) to solve the integration. The functions signature should look like this, where the final argument (‘dur’) indicates the duration of the trajectory to simulate.
def trajectory(init_n, init_y, init_c, dur) :
…
Tips:
• Experiment to pick an appropriate time step value.
• Be careful to pass time (and not the time step or the iteration count) to dpdt().
• Make your code as clear and understandable as you can to improve chances of getting partial credit.
In [ ]:
## Complete this to answer Q5
def trajectory(init_n, init_y, init_c, dur) :
…
Step 6: Plotting results and interpretation¶
Simulate a trajectory that starts with $n_0 = y_0 = p_0 = 0$ and that runs for 45 time units.
Q6a. On a labelled figure, plot a time-series that shows how the population density changes over the course of the simulation.
In [ ]:
Q6b. On a separate labelled figure, plot a time-series that shows how the temperature changes over the course of the simulation.
In [ ]:
Q6c. Superimposed on the same figure as Q6b, plot a time-series that shows how the temperature would have changed if the daisies had no influence upon temperature. Label which line is which using a legend.
In [ ]:
Q7. Interpret the results. What has happened? What are the striking features about the trajectories? What affects them? What conclusion(s) can we draw from the simulation? Write a short paragraph summarizing your interpretation.
In [ ]: