Hi, modelers.
I've been really enjoying the discussions around module two.
The tools we're seeing right now are completely different
from those we saw in module one.
We are looking at cellular automatons
that have a lot of structure and also random walks,
which are now stochastic processes,
meaning that when you rerun them,
you can get different outcomes,
because you're flipping random coins
and making random decisions.
That's very different from compartmental models
that had barely any structure in them,
or ordinary differential equations
that were deterministic by nature.
You could rerun them and always,
and always get the same results.
And so what I'd like to do in this video
is show you that even though those two different toolboxes
may seem very different, they can sometimes join forces
and help you answer complex questions.
So the project we'll be looking at
is "How the Forest Got Its Shape."
That's one of those classic scientific questions
that appears super simple, but it's hard to answer.
It's a question that popped in my mind around 2014
when I joined the Santa Fe Institute.
I saw a seminar by ecologists
who were talking about all the great data science tools
that their field was developing.
What you see here on the screen is one of them.
So this is from a paper by Hansen
and colleagues in "Science" in 2013,
and they show these high resolution global maps
that not only give us great spatial information
about forests, but also temporal information.
So the four panels here are four different places
in the world, so regions of Paraguay in panel A,
Indonesia in panel B, United States in C, and Russia in D,
and they not only show you that there's a great diversity
in terms of forest shape and structure,
so the forest would be the really bright green pixels here,
but also different dynamics across the globe.
So you also see forest loss in the red areas,
forest gains in the blue areas,
and a mixture of both in purple.
And so during the seminar, I was looking at the diversity
of forest shapes and patches and different dynamics,
and I was trying to make sense of it.
So, I raised my hand and I say it a lot in class,
there are no dumb questions, and I just asked,
"Well, what is the shape of a typical forest?
"Is what we're seeing, the ridges or rivers we're seeing
"in (indistinct) that typical
"to have linear forests like that?"
Or, "Are most forests like what we're seeing in panel B
"and have this dense shape, almost like discs or squares,
"which you would expect the forest to grow outward?"
And really, the ecologist didn't know,
and it sparked a great collaboration.
And the questions we were trying to answer
is not only what is the typical shape of a forest,
'cause you can do that by just looking
at these high resolution maps.
You can do an observational study on the data itself.
But what we wanted to answer
is really how the forest got its shape,
and this how question gets at mechanisms.
What are the key mechanism
behind forest growth and forest loss?
And if we understand these mechanisms,
maybe we can better understand forest stability
and contribute to the strength of forest ecosystems.
And so to answer that question, we're going to focus
on a subset of these high resolution maps.
We're gonna look at the Brazilian Cerrado.
It's interesting and an opportunity for two reason.
One is that it's just a very important region.
So the Cerrado borders the Amazon rainforest,
which, as you may know, is the long of the earth.
But it also has this really interesting feature
that it's a mix, so what we call an ecotone.
When I say we, I learned that word during that project.
An ecotone is an interface
between two different ecosystems or environments,
in this case between rainforest and grassland,
like savanna, like really dry grassland.
And that's interesting to see these two different systems
cohabiting under the same roof in the same region,
because it makes you wonder whether the Amazon rainforest
could be destabilized and collapse to savanna,
which is just much more dry
and obviously, less of a lung for the planet,
so that's an important question to answer.
But also, one interesting thing about the Brazilian Cerrado
are all the protected areas shown in red
where there is limited human activity,
and so we're studying forests
in their natural environment, if you will,
with not too much human interference.
So this is almost a roadmap for what we wanna do.
So we want to go from reality which is here,
represented by a picture of the system we care about,
and here, you see this mix of trees
but also very, very dry grass.
And we want to go from this reality to the data,
which is what's shown in panel B.
So this is one of those high resolution maps of tree cover
where every pixel in this case is about 30 meters
by 30 meters, so not quite like individual trees,
but enough to get a high resolution map
that gives us the shape of the forest patches.
And this is inferred tree cover,
which is a continuous variable between zero and 100%,
telling us how much of the canopy do we think
is closed by the trees.
But in most cases, it's actually pretty bimodal.
It's either green, so high tree cover
or yellow, very low tree cover,
so it feels natural to binarize it,
and that's where we already get close to a CA
at least in terms of image in panel C.
So this panel C is a finalized version of panel B
where we just apply a threshold.
Everything that has more than say 50% tree cover,
we call forest, everything that's below, we call grass,
and then we get this binarized version
of the data in panel B.
And then what we want to do is build a model
that will grow synthetic forest,
and that's what's represented in panel D.
And hopefully, the statistics are similar between C and D,
such that we hope we've unlocked some of the key mechanisms
that govern forest growth in this area and therefore,
forest stability to help us understand the stability
of the Amazon rainforest.
All right, so we wanna build our system.
We're gonna follow our recipe,
so the first big picture questions about the system itself.
What is our system?
In this case, we care about the forest grassland,
ecotone, the interface of forest and grass,
so we're gonna ignore water, like river and so on.
We're gonna ignore topology, like hills.
We're gonna ignore humans. We're in protected areas.
But we're also gonna ignore animals,
which can play a big role in shaping the growth of forests
through seed dispersion or destruction of young trees.
Is the system open or closed?
So we're gonna focus on a grid, like in the previous image.
So we're gonna assume that we're gonna ignore everything
that happens outside of that grid, so it's a closed system.
The topology or the space of the system, again,
is enforced by the data that we have.
We wanna build a system and then compare it to the data,
so we're gonna use the same topology,
which is a square grid of pixels,
30 meter by 30 meter for each pixel.
Is the state of the system discrete or continuous?
It's gonna be discrete.
We're gonna employ a cellular automaton,
so discrete state for every cell in the grid.
Is time discrete or continuous?
This is a tricky one, and here, I'm gonna answer both,
and I'm gonna add a question
that I rarely ask in those terms.
What are the important timescales,
because that's why I have this weird answer
to question number five.
So there's a tricky thing here,
which is that we have different processes
that occur over a very different timescale.
A forest fire burns in a matter of hours or days hopefully,
whereas tree growth usually happens over decades.
And so, if you think about the type of temporal resolution,
either in a step of an integrator for an ODE
or time windows for a discrete model
needed to fall forest fires,
it's gonna take that discretization,
a very long time to follow tree growth.
If you're looking at minutes per minutes
and following something that happens over a decade,
it's just gonna be a very slow model.
So that's why we're going to be combining
different temporal scale in the same model
and different notions of time.
And so now we have our basic setup,
and we want to think about our model in a compartmental way.
So even if we're building a CA,
it's useful to think of our model in terms of parts.
So boxes, what are the different types of elements
in our system or the different states
that the cells in our cellular automaton can take?
In this case, we're gonna use four.
So it might be useful for you all
to visualize the diagram in terms of boxes and arrows.
As I present this, we're gonna have four boxes.
So cells can either be occupied by trees,
so forest cell, or grass or fire or ashes.
And the next step is obviously to draw mechanisms
or arrow between all of these states.
So we'll start with the final state, ashes.
So ashes is what comes from fire,
and what can happen to ash sites?
So in our cellular automaton,
we're gonna have some sites that are just ashes
that are assigned a specific value,
and then what kind of transition can they undergo?
We're gonna say that ash sites can turn to grass
spontaneously at a rate lambda,
so that's just grass being reborn
from the ashes of the forest
or the grassland that was there before.
But it can also turn to trees
by spreading mechanism at rate alpha.
And when I say spreading mechanism for trees
and not for grass, what I mean is that
regardless of whether there is grass around or not,
grass can come from ash,
but trees will only spread to ash sites
if there is a neighboring tree.
So essentially for every ash site,
you would look in your Moore or von Neumann neighborhood,
and you would ask,
"Is there a tree there in this neighboring cell?"
If yes, there's a certain probability,
alpha times time window delta T,
that my ash site turns to tree.
What happens to grass sites?
Well, they can catch fire hereby by lightning strike,
so which will happen at a rate f or by spreading.
So if there is fire in a neighboring site,
it can spread to the grass,
and we're gonna use a probability p subscript capital G,
meaning the probability of fire spreading in grass.
If the grass sites don't catch fire, it can turn to trees,
again, either by seeding at rate alpha,
so the same mechanism as before,
but here also through spontaneous seeding at a rate beta.
This mechanism is meant to represent a bird carrying seeds
over a long distance.
This mechanism doesn't happen in ash sites.
This is all expert advice.
The idea being that to protect the seed
and help foster the growth of the tree,
it's good to already have some vegetation present.
So a lot of this is not just me
coming up with a model out of nothing
but is discussions with experts
to decide what mechanisms are at play here.
The tree sides are what's really interesting.
So now we have a way to catch fire, spread fire,
regenerate grass and all that.
The only thing that trees will do
is catch fire by spreading.
So a lightning strike won't light up a tree.
The idea is that this is rainforest territory,
so there's a lot of humidity under the tree,
that their canopy protects the humidity
or traps the humidity.
So the only way the trees can really catch fire
is if there's a fire in a neighboring cell,
and then it's spread with some probability p
under subscript capital T.
And then we assume that this PT
is much smaller than this pG,
'cause dry grassland catches fire
much more easily than trees.
And then fire sites will turn to ash with probability 1,
just meaning that fire lives for one time unit
in the cellular automaton.
So I call it the cellular automaton.
Again, it's a stochastic one, but it is a fixed set of rules
and a finite number of states that live on a discrete grid,
and that gives us a simple way to grow trees.
So if we apply it with some parameters informed
by our friends and experts in forest ecology,
we might get something like panel D here.
And here, I've drawn everything that is not a tree,
so whether it's ash or or grassland in yellow
and then trees in green to again, represent the binarization
that we did over the data.
The data itself doesn't tell us
whether it's actually grass or not.
We just assume that everything that was not trees was grass.
Then we're gonna look at a few observables.
We now have the rules.
If we think about our recipes in terms of IOU,
going backwards, we've defined the rules
that allow us to update the system.
The observables that we're gonna be looking at
is the total fraction of cells that are covered in trees,
but we're also gonna look at individual forest patches,
meaning patches of neighboring forest cells
that are interconnected, and then we're gonna look
at the area and the parameter of these forest patches.
So this is our O in the IOU, what we observe.
And typically for initialization,
we're gonna start with random initial values
and let the forest evolve on its own.
So we start with random initial values.
We start with parameters informed by the literature,
but not fine tuned, and here's what we get
in terms of comparison with the data.
So on the left what you have is a scatter plot
of the area on the horizontal axis
and meters squared for every forest patch
and their parameter in meter.
I apologize, the y label here is wrong,
parameter in meters for the same forest patch.
You have six sets of data points on the scattered plot.
So in circles, you have three subsets of the data,
depending on whether they're in the bottom third
in terms of yearly precipitation,
medium third or highest third.
And you have three sets of simulation
with somewhat different parameters,
but none of them were particularly selected
or fine tuned, I should say.
Well what's crazy about this plot
is that there is very little scatter.
Most of our simulations and most of the actual data
fall very steadily along the same axis,
meaning that there is, because we're in a log plot here,
there's a power log relationship
between the perimeter of a forest patch and its area.
And it's interesting to see that the shape of the forest
doesn't seem to depend on our parameters,
the same way that it doesn't seem
to depend too much on precipitation either.
Really, what might be important
are the mechanisms underlying the growth of these forests,
not the exact parameters.
On the right, we show the probability distribution function
for one of these dimensions, so for forest patch area.
And here again, you see a parallel distribution
or roughly parallel distribution,
because we have this very straight decay
of forest patch area as they grow bigger.
And again, very similar values
regardless of precipitation value or simulation.
I don't want to go too much into the detail here,
but the crazy thing is the fact
that we get this similar outcome
regardless of the parameters we put in.
Then we can look at our global observable,
the total fraction of cells covered in trees.
I showed two different set of parameters here.
So in both panels on the horizontal axis,
what I'm varying is the probability of fire,
which was this f rate in my rules.
So that's how often grass sites catch fire.
Notice that we have very different scales here.
On the left, it's very rare fire strikes,
and on the right, is very frequent fire strikes.
They differ by about 10 to minus six
or a factor of a million.
And then on the horizontal axis,
we look at the steady state value of tree cover,
so what fraction of the system is covered in trees.
If we start with the right panel, panel B,
then we know what to expect if there's no fire at all,
because I said the only thing that can happen to trees
is to catch fire if a neighboring cell is currently on fire.
Well, that means that if trees only have trees around them,
they're not ever gonna change state,
and if the entire system is covered in trees,
nothing's gonna happen, that's a fixed point.
So I know that for Pf equals zero,
I have a fixed point at historical zero.
And then as I start having more and more fire
during the evolution of my system,
eventually I destabilize my fully-forested state,
T equal 1 and then my steady state or equilibrium value
of tree coverage starts falling,
and of course, that still varies in time.
So here I've colored one of the data point in green,
and that's what you see in the inset
is the time series for that data point.
So you still have fluctuation,
because it is a stochastic process.
You still have fluctuations of the tree cover,
the fraction of cells covered in trees,
but there is steady values, in this case around 85% or so.
If you go on the left, you get a very different picture.
The one on the right, by the way,
like looks a little bit like epidemic transitions
that we saw in the SSIS model or SIR model
if we were looking at the steady state value
of susceptible individuals.
If the disease or the fire is very weak,
everyone's susceptible, everyone's happy.
If the disease is strong enough, it starts spreading,
and you get this bifurcation.
On the left, we have a very different picture.
So we still have that for very weak fires,
so Pf equals zero,
you expect a steady state at a fully forested system.
But what happens here is that it actually seems to decay.
And when I say decay,
I mean that the forests get destabilized
as Pf increases in a discontinuous way.
So around Pf equals four times 10 to minus six,
you get this sudden jump,
where instead of expecting a fully forested state,
you might expect only a few percent
of your system to be covered in trees,
and that's already very different.
And there's a weird shading effect on my markers here,
because what the shading represents, the alpha value,
the transparency of my markers represent the probability
that I end up either on top or at the bottom.
So if my marker is very weak,
that means there's a very low probability
that my system settles at a steady state
of about 3% or 5% of system size.
And I'm talking about roughly Pf equals four times 10
to minus six and there is therefore very high probability
of ending up in the fully forested state.
Pf increases the probability switch
and now I'm much more likely to be in the lower branch
by the time I reach an equilibrium.
What's interesting is that
if we look at a typical time series again in the inset,
we now see very large fluctuations,
so fluctuations of orders of magnitude
more than the actual equilibrium value.
And that's because sometimes, stochastically, forests grow
until there's a fire big enough to shrink them back again.
So we wanna understand this.
It's very hard to study bifurcation diagrams like this
and stochastic processes
and systems within an absorbing state.
So as soon as our system is fully covered in trees,
it's frozen there forever,
it's hard to study the fluctuations,
just a lot of computational firepower is needed.
So we can turn back and say,
"Well, what if we had built differential equations
"from our boxes and arrows
"this same way we did in module one?"
Of course now, we're not gonna be able
to describe the shape of forests, but we might be able
to better understand the stability of our fixed points.
So instead of looking at stochastic dynamics,
we're gonna fall back on mean field descriptions.
Again, mean field means here we're gonna ignore space
and simply look or assume
that everything is potentially connected,
meaning that a cell is defined by its state
but then connected to the average
or expected value of all neighbors.
So let me describe just one question in detail,
and here, I'm switching some of the rates a little bit,
because it's almost impossible to get a one-to-one matching
between parameters in a stochastic and discrete simulation
and parameters in a mean field description,
so some of them are gonna be a little different,
and you'll see.
So I'm just gonna describe this first differential equation,
which is B dot, so the time derivative
for the number of burning trees.
So the number of burning sites,
so that could be burning trees or burning grass,
so burning sites in my lattice, in my grid.
So every grass site, G,
can catch fire with probability f, or with rate f.
So I'm gonna increase of f times G
in the number of burning sites.
That's one arrow flowing from G to B.
I also have another arrow flowing from G to B,
which is the probability that grass sites
have a neighboring site that is currently burning
and that the fire spread.
Why is there a factor of four here?
Well, let's take the perspective of a single grass site.
So for every grass site here, we have four neighbors
that we're specifying our neighborhood here.
So we have four neighbors,
which will be burning with probability B,
because there is no structure,
so every neighbor is a random neighbor.
And for every of those potential
four times B burning neighbors,
the fire will spread with probability row underscore G.
And here we use row, because it's not a probability,
I just misspoke.
It's not a p underscore G as in the cellular automaton.
This is an actual rate, so we're using row,
because it is a rate per time unit, not a probability.
Then we have a very similar terms for trees catching fire.
So for every tree site T, we have four neighbors
that are currently burning with probability B,
therefore an average four times B burning neighbors,
which spread the fire with rate row subscript, capital T.
And then the only arrow leaving the burning box
is the arrow leading from burning B to ash A,
and that happens with a rate here that we're calling mu,
which was previously one,
but here we have to put an actual value on it.
It's just a very fast rate
at which local fires stop burning.
And you can go through the same process
and understand every term in the growth of the system.
You'll notice that here in the grass term,
we're using exponential growth and not logistic growth.
That's a choice, assuming that at any given time,
there won't be a lot of ash present in the system.
And then, we can solve
for the bifurcation diagram of the system.
I'll focus on panel A for now.
So we have a lot of parameters, so it's not as easy
to analyze as the SSIS model, for example.
We're gonna be looking at panel A at the 3-D plot,
so the steady state or the fixed value of forest coverage
as a function of fire rate on the x axis.
On the depth axis, so let's call it y,
we're gonna be looking at the growth rate of trees beta,
which if you remember, is the seeding mechanism for trees.
So this is how much do birds
from outside of our system carry seeds inside?
And that's why I said that we sort of capture an open system
but not really, 'cause we don't model the birds
or what's out there, we just have this parameter
that accounts for that.
And so what we see here,
the blue lines are the usual fixed points.
So, we looked at that in terms of SSIS before.
It's the the expected fixed point
for a given value of f and beta.
What happens in red
is that sometimes those fixed points are unstable.
So what connects the fully forested state T equal one
to the lower branch of stable values
is a branch of unstable fixed points.
We've seen unstable fixed points
when looking at flow diagram before,
so that's for a given set of parameters.
Here, we're changing the parameters along the x and Y axis,
and therefore, our fixed point is moving
and tracing this red line.
When it becomes wide here,
it's simply because the red line goes into a physical value
for fire rate negative and stuff like that,
so it's just not a region of parameter space
that we wanna care about, so we're just not plotting it.
But, they're still connected, and they're still out there.
Well, what's interesting here
is that if we have a lot of spontaneous tree growth,
which again, could be birds,
it could also be me going out and planting trees,
we're stabilizing the system.
We're helping not only the lower branch
of stable values reach higher levels of forest coverage,
but we're getting rid of this discontinuous instability,
and that's good news,
because these discontinuous instability
are what we're really scared about,
is for the Amazon rainforest
to suddenly go from a highly forested state
to a barely forested state, a sudden discontinuous collapse.
What we want is a stabilized system.
So one key mechanism already that we could suggest
is spontaneous tree growth,
either by replanting in grassland
or by increasing biodiversity in the system.
All of these would be a good option.
We can also look at different ways in terms of network.
We can expect a disconnected transition
or continuous transition, some cases, no transition at all,
meaning that if trees grow quickly enough,
there's no fire strong enough to destabilize the forest
and then we can solve for all that.
I won't go too much into detail.
I'll give you the reference,
but this is a good example of chapter eight on bifurcation.
There's a lot of rich dynamics here at play.
The reason I wanted to talk about this project
is that we can also make interesting predictions.
Once we believe that we've captured something
about the growth of forests in this system,
we can make predictions about individual forests.
So this is a cool experiment that we did.
So we take our realistic set of parameters,
we run our cellular automaton,
and at some point we freeze the system,
and then we go in the system with scalpel,
and we find a patch of grassland somewhere,
and in that patch of grassland,
we inject a forest with a very specific shape.
So we try all possible shapes,
meaning that on the x axis we have area
in terms of number of cells.
So if you have a forest that covers an area of 100 cells,
it could be a line with a parameter of 100 on one side,
100 on the other, and two on both ends,
so a parameter of 202, or it could be a square
that's 10 cells by 10 cells,
you can get very different shapes,
and we're gonna try them all essentially
over this little figure.
The black lines in these figures
are showing what shapes are possible.
So you can't be any more compact,
meaning you can't have a parameter lower
than what you would have as a square.
So that gives us this square root relationship
between the perimeter of the forest patch and its area.
And you can't have a perimeter higher
than if you were a line.
And so that gives us this linear relationship,
which is the top black line that is steeper,
so the top black line as a slope of one,
the bottom one is slope of one half, a square root,
and all possible shapes exist
in between these two lines, so in this cone.
In the red, you have the stable prediction
that we make from our model,
which is an exponent of about three quarters,
so meaning a perimeter
that's about three quarters the area of the forest patch.
Then we try a bunch of individual patches,
we see how they grow, how does their area grow,
how does their perimeter grow,
and we average it over thousands and thousands
of realization of the same forest patch.
So here, are probably billions and billions
of simulations of individual forest patch
within a larger ecosystem.
And the orange arrows are showing you the expected growth
or expected dynamics,
because sometimes it's not growth at all,
the expected dynamics of the forest patch.
So if you start below the red line,
you get forest patches that grow a tiny bit in area.
But remember, growing an area is a slow process
and they grow a lot in perimeter.
Growing in perimeter means that you're facing fire,
and the fire is chipping at you,
so you go from a square to something with a (indistinct)
and a higher perimeter, so that's what happens.
And as they do so, the forests grow in perimeter,
but the arrows get smaller and smaller.
So around the red line, there's actually an area
where forests barely follow any dynamics.
They're very stable. They grow as much as they get burned.
But if they burn too much,
which will always happen at some point,
because it's a stochastic process, if they do burn too much,
then they enter in this, I like to call it, it's like river,
but this flow of rapid decrease in both area and perimeter
that happens at the top.
So once you're too linear, what happens,
that fire can fragment forest patches
and just make them collapse really, really quickly.
So stochastically, you'd expect forests
to slowly go towards the right
and up following the red line,
but if they're unlucky and there's a lot of fire,
then they can suddenly collapse.
That's the lifecycle of a forest patch.
What's cool about this is that
if you were to go cut down some trees,
deforestation, part of an industry,
well at least this gives you a map
of what shape you should leave the final patch in.
Or if you're tempted to leave the forest patch
in a weird shape because it's convenient,
this can give you an idea of the risk
that you're putting this forest patch in.
Here's a reference, the paper I wanna point out.
This was a great multidisciplinary collaboration
between physicists and ecologists and more.
I want to point out, especially the team
that was with me at the Santa Fe Institute,
so Andrew M. Berdahl, we're gonna hear more about him
who does a lot of collective behavior ecology.
Sidney Redner, I already recommended Sid's book,
on "A Guide to First Passages Processes,"
one classic if you're interested in random walks
and processes like this.
Uttam Bhat, which was with us at the Santa Fe Institute
and then Adam Pellegrini and Stephen Pacala,
who were at Princeton at the time in ecology.
So, there are a lot of directions this could take,
like once you have a good model,
there are a lot of cool questions to ask.
We could think about comparing data
within all those protected areas that we looked at
and outside to see what is the impact
on the shape of forests of human-