Ecological Dynamic Programming Optimization in Python

It’s been awhile since I’ve posted anything. I’ll use the excuse that I’ve been busy, but mostly I just forget. Regardless, I’m learning how to do Dynamic Programming Optimization (DPO), which sounds more complex than it is. The reason for this is that DPO allows us to simulate the behavior of individuals who make decisions based on current patch quality.  DPO is an exciting tool that forms the foundation of individual-based models, which allow us to assess community and ecosystem dynamics as emergent properties of individual decisions based on well-grounded, physiological principles (hence my interest).The underlying idea, as I understand it, is that individuals assess their future reproductive output prior to making a decision (I’m interested in optimal foraging, so I’ll put this in the context of foraging). We can model how individuals make decisions over the course of a lifetime, and track each individual, which then allows us to make quantitative statements about populations, communities, or ecosystems.

This all sounds complicated, and it can be difficult. So it’s easiest to jump right in with an example. Here’s the code, along with a basic explanation of what’s going on, for the very first toy example found in “Dynamic State Variable Models in Ecology” by Clark and Mangel. In the book, they give the computer code in TRUEBASIC, but… in all honesty.. no one I know uses that. I could use R, but we all know my feelings on that. So here’s their toy model programmed and annotated in Python.


Suppose we’re interested in the behavior of fish and how they choose to distribute themselves among three different patches. Two patches are foraging patches and one is a reproduction patch. The first thing we need to do is make an objective fitness function F. In this example, F relates body size to reproductive output. It can be thought of as what the organism believes about its final fitness given a particular body size. Let’s make it increasing, but asymptotic, with body mass:

F = \frac{A(x-x_c)}{x-x_c+x_0}

Here, A is sort of a saturation constant beyond which fitness increases minimally with body size, x_c is the body size at which mortality occurs (0), and x_0 is the initial body size. Lets set x_c=0, A=60 and x_0=0.25*x_{max}, where x_{max}=30 is the individuals maximum possible body size. You have to set the ceiling on body size otherwise organisms grow without bounds.

OK so that’s what the organisms “believes” about its fitness at the end of a season given a specific body mass. The question is: How does the organism forage optimally to maximize fitness at the end of the its lifetime?

The obvious answer would be to simulate a foraging decision at every time step moving forward, and then have it decide again at the new time step, etc. This is computationally expensive, so to circumvent this we work backwards from the end of the season. This saves time because there are far fewer paths to a known destination than there are to an unknown one (essentially).

So we imagine that, at the final time step t_f, the organism’s fitness is given by F for any given body mass.

import numpy as np
import matplotlib.pyplot as pet
import pandas as pd

n_bodysize = 31 # 31 body sizes, from 0-30
n_time = 20 # 20 time steps, 1-20
n_patch = 3, # 3 foraging patches

# make a container for body size (0-30) and time (1-20)
F = np.zeros((n_bodysize, n_time))

# make the function F
x_crit = 0
x_max = 30
acap = 60
x_0 = 0.25*x_max 
t_max = 20

# calculate organism fitness at the final time point
for x in range(x_crit, x_max+1):
    F[x,-1] = acap*(x-x_crit)/(x-x_crit+x_0)

Now that we know fitness at the final time point, we can iterate through backwards (called backwards simulation) to decide what the optimal strategy is to achieve each body mass. To determine that, we need to know the fitness in each patch. Let’s start with the two foraging patches, Patch 1 and Patch 2. We need to know four things about each patch: (1) the mortality rate in each patch (m), (2) the probability of finding food in each patch (p), (3) the metabolic cost of visiting a patch (a), and (4) the gain in body mass if food is successfully found (y). For these two patches, let:

Patch 1: m=0.01, p=0.2, a=1, y=2

Patch 2: m=0.2, p=0.5, a=1, y=4

Right away we can see that Patch 2 is high-risk, high reward compared to Patch 1. In each patch, we calculate the next body size given that an animal does (x‘) or does not (x”) find food:

x' = x-a_i+y_i

x'' = x-a_i

Those are simple equations. Body size is the current body size minus the metabolic cost of foraging in the patch and, if successful, the energy gain from foraging. Great. Now we can calculate the expected fitness of each patch as the weighted average of F‘ and F” given the probability of finding food, all times the probability of actually surviving. For these two patches, we make a fitness function (V):

V_i = (1 - m_i)*[p_i*F_{t+1}(x') + (1-p_i)*F'_{t+1}(x'')]

The reproductive patch is different. There is no foraging that occurs in the reproductive patch. Instead, if the organism is above a critical mass x_{rep}, then it devotes all excess energy to reproduction to a limit  (c=4). If the organism is below the reproductive threshold and still visits the foraging patch, it simply loses mass (unless it dies).

OK this is all kind of complicated, so let’s step through it. We know what fitness is at the final time step because of F. So let’s step back one time step. At this penultimate time step, we go through every body mass and calculate fitness for each patch. Let’s do an example. If x=15, then we need to know fitness in Patch 1, Patch 2, and Patch 3. For Patches 1 and 2, we need to know the weight gain if successful and the weight gain if unsuccessful.

x' = max(15-a_1+y_1, x_{max})

x' = 15-1+2 = 16

x''=min(15-a_1, x_{c})

x'' = 15-1 = 14

The min and max functions here just make sure our organism doesn’t grow without limit and dies if metabolic cost exceeds body mass. So these are now the two potential outcomes of foraging in Patch 1 given x=15. The expected fitness of these two body masses is given as F(16) and F(14). Plug all these values into equation V for Patch 1 to get the expected fitness of Patch 1 at a body size of 15. We then take the maximum for all three Patches, save whichever Patch corresponds to that as the optimal foraging decision, and then save as the fitness for body size 15 at that time step. So at body size 15, for Patch 1 is  39, for Patch 2 is 38, and for Patch 3 is 37.6, so the individual will foraging in Patch 1 and now fitness for body size 15 at this time step is 39.

Repeat this procedure for every possible body size, and you’ll get the fitness for every body size at the second to last time step as well as the optimal foraging patch for every body size at that time.

Then, step backwards in time. Repeat this whole procedure, except now the value for each body mass doesn’t come from the equation F, but comes from the fitness we just calculated for each body mass. So for example, if an organism’s foraging decisions at this time step lead it to a body mass of 15, then F is now 39. Again, repeat this for every body mass, and then step back, etc.

Here’s the full Python code for how this is done:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

n_bodysize = 31 # 31 body sizes, 0 - 30
n_time = 20 # 20 time steps, 1 - 20
n_patch = 3 # number of patches

# make a container for body size (0-30) and time (1-20)
F = np.zeros((n_bodysize, n_time))
# make a container for three patches, each with body size (0-30) and time (1-20)
Vi = np.zeros((n_patch, n_bodysize, n_time))
# make a container for the optimal decision
d = np.zeros((n_bodysize, n_time))

# make a container for the mortality rates in each patch
m = np.zeros(n_patch)
# make a container for the probability of finding food in each path
p = np.zeros(n_patch)
# make a container for the metabolic cost
a = np.zeros(n_patch)
# make a container for food gain in each patch
y = np.zeros(n_patch)
# make a container for reproductive output in each patch
c = np.zeros(n_patch)

x_crit = 0
x_max = 30
t_max = 20
x_rep = 4

m[0] = 0.01; m[1] = 0.05; m[2] = 0.02
p[0] = 0.2; p[1] = 0.5; p[2] = 0
a[0] = 1; a[1] = 1; a[2] = 1
y[0] = 2; y[1] = 4; y[2] = 0
c[0] = 0; c[1] = 0; c[2] = 4

acap = 60
x_0 = 0.25*x_max

# Calculate Fitness for every body mass at the final time step
for x in range(x_crit, x_max+1):
    F[x,-1] = acap*(x-x_crit)/(x-x_crit+x_0) # maximum reproductive output for each body size at the final time

for t in range(0, t_max-1)[::-1]: # for every time step, working backward in time #print t
    for x in range(x_crit+1, x_max+1): # iterate over every body size # print x
        for patch in range(0, 3): # for every patch
            if patch in [0,1]: # if in a foraging patch
                xp = x-a[patch]+y[patch] # the updated body size
                 xp = min(xp, x_max) # constraint on max size
                 xpp = int(x-a[patch]) # updated body size if no food
                 Vi[patch,x,t] = (1 - m[patch])*(p[patch]*F[int(xp),t+1] + (1-p[patch])*F[xpp, t+1]) # Calculate expected fitness for foraging in that patch
                 if x < x_rep: # in a reproduction patch
                     xp = max(x-a[patch], x_crit)
                     Vi[patch, x, t] = (1-m[patch])*F[int(xp), t+1]
                     fitness_increment = min(x-x_rep, c[patch]) # resources devoted to reproduction
                     xp = max(x-a[patch]-fitness_increment, x_crit) # new body size is body size minus metabolism minus reproduction
                     Vi[patch, x, t] = fitness_increment + (1-m[patch])*F[int(xp),t+1]
        vmax = max(Vi[0,x,t], Vi[1,x,t])
        vmax = max(vmax, Vi[2,x,t])
        if vmax==Vi[0,x,t]:
                 d[x,t] = 1
        elif vmax==Vi[1,x,t]:
                 d[x,t] = 2
        elif vmax==Vi[2,x,t]:
                 d[x,t] = 3
        F[x,t] = max # the expected fitness at this time step for this body mass

This model doesn’t really track individual behavior. What it does is provides an optimal decision for every time and every body mass. So we know what, say, an individual should do if it finds itself small towards the end of its life, or if it finds itself large at the beginning:

T,X = np.meshgrid(range(1, t_max+1), range(x_crit, x_max+1))
df = pd.DataFrame({'X': X.ravel(), 'T': T.ravel(), 'd': d.ravel()})
df = df.pivot('X', 'T', 'd')
sns.heatmap(df, vmin=1, vmax=3, annot=True, cbar_kws={'ticks': [1,2,3]})


And that’s that!! To be honest, I still haven’t wrapped my head around this fully. I wrote this blog post in part to make myself think harder about what was going on, rather than just regurgitating code from the book.

Saturation Science

The other day I introduced our upcoming Aquarius mission, so I thought today I would take a few minutes to introduce the science behind it all. While the opportunity to live underwater and dive nine hours a day seems like reason enough to me, it can get pricey to pull off such an endeavor and “because it’s awesome” typically doesn’t get a project funded. In my previous post I mentioned  I’d be studying herbivory and grazing across a depth gradient so I thought I would briefly discuss why it’s important to understand.

I won’t spend much time here arguing why reefs are important; suffice to say, if you don’t appreciate the intrinsic value of a reef, or that ~ 10% of the world’s population depends on them for food, you may appreciate the services they provide, such as tourism, fisheries, coastal protection and natural products, which could be worth as much as $375 billion per year. Despite their various values, coral reefs are in bad shape these days and their prognosis looks bleak. Warming oceans, overfishing, changes in seawater chemistry, and more severe tropical storms are just a few of the factors that have contributed to the loss of as much as a third of world’s coral reefs in just a decade. When these slow-growing corals die on shallow reefs, faster- growing algae quickly take their place. These algae smother and kill the remaining corals, while preventing new coral from settling. The job of preventing this transition from a coral- to algae-covered reef falls to coral reef herbivores such as urchins and parrotfishes, which eat algae and maintain the open space for swimming coral larvae to settle and grow. Unfortunately, due to human impacts the speed at which corals are now dying is so fast that herbivore populations, already decimated by overfishing and disease, are often unable to keep pace with the algae growth.

A time series taken by geologist Eugene Shinn showing the same location in the Florida Keys from 1959 - 1998.

A time series taken by geologist Eugene Shinn showing the same location in the Florida Keys from 1959 – 1998.

But there’s still hope for the future. In addition to traditional strategies, such as better management of shallow reefs, deep-water or “mesophotic” coral reefs that grow between 30 m and the depths that light can no longer penetrate, are gaining attention from researchers as naturally protected “source” populations for surrounding areas. Because the deep waters above these reefs can provide a natural buffer from threats such as artisanal fishing, storm damage, and the high-temperatures that trigger coral bleaching, these reefs are naturally protected from many of threats humans have created for shallower reefs. As a result, these deep reefs may act as a sort of self-sufficient Marine Protected Area, where species survive and their offspring eventually repopulate the surrounding shallower reefs. Due to the technical difficulties of working at depth, researchers are only recently beginning to understand the role of these deep coral reefs and the processes important to their continued health. While we know herbivory is important for shallower reefs, declining numbers of herbivores at depth, coupled with slower rates of algae growth, have created uncertainty in how important herbivory is to maintaining the health of deep-water coral reefs. And that’s where our science fits in. Starting on Sunday we will begin documenting herbivore communities, feeding patterns and algae growth to understand the role of herbivores in structuring these important systems.

The future of the forest in the clouds

Increasing temperatures associated with global climate change have resulted in subtle but noticeable changes in the way species act in both time and space. For example, spring flowering and budset for many plants has come earlier and earlier, birds migrate and nest earlier in the year, and insect outbreaks have become more severe. One area where climate change is predicted to have an especially dramatic impact is in tropical mountain cloudforests.

Unless you’ve spent a lot of time trekking around the Andes or retracing Alfred Wallace’s steps in Indonesia you have probably not noticed that most tropical mountain cloudforest trees are already responding to climate change. These ecosystems are highly diverse, unique ecosystems because the forest is constantly engulfed in clouds and fog, creating an eerie yet fascinating ecosystem where life is literally clinging to every surface.


Inside a tropical mountain cloudforest where every tree trunk and branch is covered in moss, orchids, and bromeliads.

On the eastern slope of the Andes this cloudforest is created when warm moist air moves westward across the lowland Amazonian jungle and runs into the very steep Andes. When this air hits the mountains it is forced upwards, cools, condenses and drops most of the moisture as precipitation or fog. Similarly, if you were to walk up this same mountain you would notice that temperatures decrease rather quickly as you ascend the mountain. What you may not notice is that the composition of the species present on the mountain also changes as you climb. Species turnover is so high because most species on these mountainsides are adapted to live within a narrow belt of temperature and precipitation, which is associated with a specific elevation. Moving just 300 vertical meters in the Andes can result in a complete change in the forest, from the species of trees to the birds and even down to fungus found in the soil.

As temperatures increase these species need to migrate upslope towards cooler, higher refugia because temperatures in their current range will quickly become too hot. This could lead to several scenarios where species shift their ranges in response to increasing temperatures. However, the most likely migration scenario is that species will slowly move upslope at the high elevation edge of their range but lose habitat at the lower edge of their range (Figure part D). This would lead to range contractions and elevated extinction risks in this hyper-diverse cloudforest ecosystem. In fact range contractions are already occurring on many tropical mountains including insects on Mt. Kinabalu in Borneo, lizards and snakes in Madagascar, and birds and trees in Peru.


The majority of tropical mountain cloudforest species are predicted to move upslope with changing climate. How they migrate upslope is important for the future conservation of this ecosystem. In scenario A the species range does not move upslope but the relative abundance of individuals within the range shifts upslope. An alternative is that the entire range, both the high elevation and low elevation edges shift upslope (B). One unlikely scenario is that the species high elevation edge shifts upslope and the low elevation edge remains the same (C). The most likely and most common scenario is that the high elevation edge of a species range does not shift upslope but the low elevation range does (D). Species will experience range contractions as the species continually lose area at the low elevation part of their ranges but do not gain elevation at the high end. Figure taken from Feeley, Rehm, Machovina. 2012. Frontiers in Biogeography.

As climate change biologists we are desperately trying to figure out how species will respond to increasing global temperatures so we can make informed conservation decisions yet our understanding of many ecosystems is so basic that we do not have much confidence in our ability to predict the future. Regardless of our predictions, we already see many species in the diverse tropical cloudforest systems suffering from increased regional and global temperatures. Unfortunately the pace of contemporary climate change far exceeds the capacity of most species to adapt to their environment and will result in major biodiversity losses in cloudforests and ecosystems worldwide.

For more information on projects we are currently working on please visit:

One big reason climate skeptics don’t believe in climate change

If you follow a credible news source at all (i.e. not Fox, CNN, or MSNBC), you’ve undoubtedly been seeing all sorts of stories: “July hottest month on record, ever!“, “Massive drought grips Midwest“, “Heat records shattered across America“, “327th consecutive month with above average temperatures” (now at least 330), etc. Yet despite these facts (indeed, temperature measurements are facts and not guesses), climate skepticism is growing in the United States! Why, in the face of all these shattered temperature records, droughts, heat waves, and snowless winters, do people still argue against climate change?

I’ve pointed out before that the delivery of the climate change message has been less than perfect; fear is not the way forward. Presenting the facts in a neutral way might possibly work, but doom and gloom scenarios will not work (even if they are true).

So what is the one big reason many skeptics still deny climate change? One word: Hypocrisy. This has been bugging me for a while. My knowledge is restricted somewhat to my experience during my Master’s and PhD. programs, but ecology graduate students really do an abysmal job of taking their own advice. For example, ecologists and climate experts say that one of the most effective actions you can take to mitigate climate change is stop driving. Carpool, bike, walk, skateboard, whatever, just stop using your car so much. That’s a pretty big lifestyle change to ask people to undertake, especially when the people giving the advice are terrible at following it themselves and they should know better.

I bike or walk to work every day. Yes, I drive to the store (sometimes) or to my friend’s house, but that’s about it. I generally use my car way less than most, and that’s still considering my friends always make me DD (I’m chalking that up to carpooling, as well). Here’s a map of a few people I worked with in my Master’s program in relation to the lab (pins are not actual addresses, but close).

I really wanted to slap the two people who drove 0.17 miles to work. For your information, that was driving from the dorms across the parking lot to the lab.

At this place, most students lived well within 4 miles of work. That’s a 20 minute ride at most at a very comfortable pace. Hell, it was an island, it wasn’t exactly a long trip anywhere.

This one is worse. Here’s a map of approximate grad student housing in relation to campus at my PhD program.


Not only are they incredibly close, but 1) It’s quicker to bike due to congested streets and the park being closed to automobiles, and 2) The cluster of people that live near one another never even carpooled, they drove separately.

Here’s another map of a place I have worked. In this situation, the lab was split into two segments roughly 3/4 miles apart. Moreover, the lab has about 30-40 free bikes for transportation around campus. Guess who would sign out a truck just to drive back and forth between labs?

If you answered ‘Marine biologists doing climate change research’, you win!

If I sound angry, it’s because I am. I’m probably going to have a few people angry at me for this (as if anyone reads this anyway). People who should know better and constantly complain about “ignorant Americans driving their stupid cars” can’t be bothered to ride a bike themselves. In that case, why should anyone believe us? Skeptics probably say, “Well, if climate change is really an issue, why are these scientists still putzing a mile down the road in a car?”

“You must be the change you to see in the world”, right? Enough ‘Do as I say, not as I’ve done’. How about, “Do as I do, because we all need to work this out together”.

Also, if you haven’t picked up on it yet, 24 hour cable news networks are awful, awful sources of real information. Also, I was careful not to say names or places in the maps above. So if you do suck at not driving, your secret is safe until someone else says it.

Why Fear is the Wrong Approach to Climate Change

Climate skeptics are growing in numbers within the U.S., currently comprising about 25% of the general population, up from about 15% or so a few years ago (goes to show that no matter how small the group, being well-funded will go a long way to hijacking the political and educational system). Regardless, there is now a growing effort to convince climate change deniers that it is real, that it constitutes a large threat to natural and human health, and that something must be done about it.

The most common approach has been to inundate the public with information about the disastrous potential of climate change on all aspects of life as we know it. Central to this approach is the idea that climate skeptics are skeptics because they lack information, and the best remedy to that is education. This doesn’t appear to be working, in fact, it may be doing more harm than good. Because the consequences of climate change are usually negative, ‘education’ usually involves telling stories about how bad things will get in the future, and right-wingers have seized on this calling it, not entirely inaccurately, ‘liberal fear-mongering‘.

Oft ignored is the psychology behind climate change beliefs. In particular, we need to understand how to get skeptics on board (which plays into the much bigger psychological picture of choices, ideology, and self-identification). A new paper in Nature Climate Change presents an incredibly interesting take on the situation. Bain et al. devised a way to experimentally assess which messages get positive reactions from climate skeptics.

They found that negative messages, essentially attempting to educate the skeptics on the effects of climate change, led to skeptics entrenching in their beliefs even further. In contrast, messages that ignored science and instead focused on their self-identity, perception within the community, and economic impacts received significantly more positive responses. That is, messages that suggested taking action on climate change would lead to a friendlier, less hostile community, more positive perception of the individual by his or her peers, and promised that taking action on climate change would create jobs actually worked! Not 100% of the time, but far more than the ‘wrath of nature’ approach.

So maybe it’s time for us to switch tacks in this debate. Science is a dirty word, and presenting scientific evidence in some camps is likely to reinforce negative opinions of climate change. Maybe now its time to suggest that taking action on climate change would help individuals be better citizens, would help this country heal as a whole, would provide capital and jobs for a huge number of green start companies. Maybe taking a ‘responsible citizen’ approach, ignoring the science, would help us make progress. On a religious front, perhaps championing the notion that we are ‘stewards of the earth’, as ordained by God in Genesis, would help climate deniers realize that our responsibility is to care for the planet (look up the definition of steward and debate its meaning later). After all, negative actions by another party only reinforces one’s belief in one’s self and ideology. If we’re trying to reach across idealogical lines, perhaps its time to be a bit more positive. Trying to be positive can’t hurt.