Interaction Networks

I took a course titled, "Mathematical Ecology" in undergrad at the University of Texas. The course was taught by Dr. Eric Pianka and it was the worst classes I ever took. Well, no, not really. It was a distressing class, though. I studied like mad and got Ds on every mid-term. I didn't realize until the night before the final how I was approaching the class wrong. Most points in each midterm were assigned on the basis of questions about a single mathematical model described in the course section (and not all the other little info-bits that I had been focusing on). I stayed awake all night, studying only the mathematical models discussed all semester. I ended up getting a B as my final grade for the semester.

I never did pick up the graded exam, but my course grade told I ended up blowing the exam out of the water. Now, this alone might not have been enough to get the grade I got. I suspect the professor acknowledged my --eventual-- mastery of the semester's most important material by making my final exam score have a bigger impact in the overall grade than it was originally designed to be. Thanks very much for that, Dr. Pianka.

---

Figure 1a. Positive
interactions.

One particular lesson started with an observation of life out in the open ocean, well away from shore or reef systems at least. Out in the open ocean, there is a large amount of biomass present in tiny plankton and in large fish, but relatively little in small fish. The plankton feeds the small fish and then the small fish feed the large fish. Another way to say this is that plankton biomass has a positive influence on small fish biomass and small fish biomass has a positive influence on large fish biomass. We can represent each positive interaction with a pointed arrow. One positive influence and a second positive influence together forms an overall positive influence (Figure 1a). This can be translated into an equation as: (I₁ * I₂) = I_f.

Figure 1b. Negative
interactions.

We can think about these steps in the opposite direction. Large fish biomass has a negative influence on the amount of small fish biomass, which has a negative influence on the amount of plankton biomass. We can represent each negative interaction with a flat-arrow, going the opposite direction of the earlier pointed arrows. One negative influence and a second negative influence together also forms an overall positive influence (Figure 1b). This can be translated into an equation as: (-I₁ * -I₂) = I_f.

Figure 1c. Fishy
feedback.

Plankton biomass has a positive influence on large fish biomass and large fish biomass has a positive influence on plankton biomass. All together (Figure 1c), we have a positive feedback loop that explains why there is more biomass in plankton and large fish than there is in small fish. On average, the ecosystem only produces enough small fish to be consumed by the large fish. Any more small fish would result in the growth of the population of large fish. Any less small fish would result in a reduction in the population of large fish.

This sort of analysis won't let you know what levels the elements of the system being modeled will stabilize at, or if it will even stabilize at all, but it can quickly give you an overview of the likely overall behavior of the system.

---

Figure 2a. Otters in
the Kelp forest.

Another scenario with the same dynamics is the relationship between Kelp, Sea Urchins, and Sea Otters (Figure 2a). This relationship is discussed in some detail in a recent episode of Science Friday. Otters help maintain the kelp forest by limiting the numbers of the major herbivore (urchins) which would otherwise obliterate the kelp.

Figure 2b. Orca in
the Kelp forest.

We can add Orca, a major predator of sea otters, into the model (Figure 2b). Because Sea Otters aren't the dominant food source for Orca, I haven't added a positive arrow between Sea Otters and Orca. The model allows us to see that Orca have a negative influence on Kelp forests. Humans hunting otters for their fur would have a greater impact than the Orca because the humans are looking specifically for the Sea Otter (no matter how scarce), while the Orca will switch to other prey if Sea Otters become scarce.

A more complicated system also discussed in the episode of Science Friday describes how the introduction of a vaccine against the Rinderpest virus could have an impact on the population of Giraffes on the Serengeti. Rinderpest is a virus of cattle which kills Wildebeest. Wildebeest eat large amounts of grass, so when Wildebeest are killed off, more grass will grow. More grass means the ecosystem is more susceptible to fires, which kill Acacia seedlings. Since adult Acacia trees feed Giraffes, more fires means less food for Giraffes. Altogether, the model presented in Figure 3 shows that the introduction of the Rinderpest vaccine would result in an increased population of Giraffes. Unfortunately, this analysis wasn't done in advance of vaccine introduction. Instead it was developed in an attempt to understand why the population of Giraffes increased after the introduction.

Figure 3. Impact of Rinderpest vaccine on Serengeti Giraffes.

---

Figure 4a. Opaque-white
switching in C. albicans.

The same logic can be used to analyze networks of interactions outside ecological modeling. Figure 4a describes the genetic regulatory network responsible for white-opaque switching in Candida albicans, a sometimes pathogenic yeast commensal of the human gut. Opaque-white switching refers to a transition of the yeast between a reproductive and vegetative developmental state. The opaque/white labels refer to how colonies of the cells in each developmental state appear when growing on media in a petri dish. My Figure 4a is updated a figure from the original paper describing this interaction network so it has colors consistent with my other figures here. The relationships illustrated in the figure had been constructed by involved experiments examining interactions between each of the genes involved in the process.

The diagram was presented by a speaker at a yeast research conference I attended. The speaker then went through an involved discussion of all the experiments they did to explain why the network resulted in the final observed behavior of switching between the white and opaque developmental states. As the presentation went on and on, I found myself getting more and more irritated. I was expecting a couple slides describing the interactions, then a continuation on to the real meat of the presentation. It was quickly clear to me why this network would result in the observed switching behavior because I was analyzing it using the techniques discussed earlier in this posting. The final conclusion of the presentation was the result I understood at the beginning of the presentation. I felt very frustrated.

After the presentation, I tried to convey my frustration to my graduate advisor. Unfortunately, I didn't know how to explain how I was processing the interaction network that had been shown. It turns out that most biologists don't have the mathematical background that I had taken for granted in my own research.

Figure 4b. Simplification step 1.

I knew how the overall network behaved because I was able to quickly simplify the network. In the left half of the network (in Figure 4a), the positive interaction between WOR1 and CZF1 and the negative interaction between CZF1 and EFG1 together have an overall negative interaction. There is already a negative interaction between WOR1 and EFG1, so this side chain would act to amplify the effect. For our purposes, this means we can simplify the network by removing this side chain (Figure 4b).

Figure 4c. Simplification step 2.

The interactions between WOR1 and WOR2 represents a feedback loop which stabilizes them in an active state. The negative interaction between EFG1 and WOR2, paired with the positive feedback loop, means that the right half of the figure has an overall negative influence between EFG1 and WOR1. Again, we can simplify the network by drawing a single negative interaction on the right (Figure 4c).

Two elements in the network with negative interactions between them represent a bistable switch. If either WOR1 or EFG1 gets the upper hand, it then suppresses the other. Only one will ever be active at one time. The final element of the network which remains is the influence of the a1/α2 ratio, which acts as an outside control element on the switch. There are other elements impacting the activity of EFG1 which aren't described in this network, but they will also act as external control elements. Together, these outside influences determine the state of the switch and when it can change from one state to the other.

---

The relationship between WOR1 and EFG1 was very quickly apparent to me at the start of the presentation, but then the presentation kept going on and on in some detail trying to explain the overall behavior. A the time of the conference, my frustration then was in why the presenter wasted so much of our time, as well as in my inability to explain this concept to my graduate advisor.

Now, however, my frustration with this memory has more to do with my irritation at researchers not collaborating with mathematicians (or at least mathematically-inclined researchers) during their work. Having some basic level of literacy at math, computer programming, statistics, or any other specialized field can help you in many ways. You don't need to be an expert in these subjects because even a limited grasp of them will help you to know when it is a good idea to call in someone with the specialized knowledge that you lack.

Dr. Pianka's course is the one I have most often thought of over the years. Specific lessons I learned in his course routinely come to mind and help guide how I interpret material seen in my scientific career so far. I expect the course will continue to inform my future biology research endeavors. Thank you again, Dr. Pianka.


References

• en.wikipedia.org/wiki/Eric_Pianka
• www.sciencefriday.com/episodes/march-18-2016/
• Kelp ecosystem:
• cmast.ncsu.edu/cmast-sites/synergy/seaotter/sskelp.html
• phys.org/news/2014-11-kelp-beds-barren-grounds.html
• en.wikipedia.org/wiki/Prey_switching
• Candida albicans:
• en.wikipedia.org/wiki/Candida_albicans
• www.ncbi.nlm.nih.gov/pubmed/17880264 

Biological Patterns: Turing and Young

The development of living creatures is all about he formation of patterns. Where hairs/feathers/scales grow on the skin; where ribs form; how many limbs grow; what colors develop and how they're organized...  Without patterns in development, none of those familiar details of organisms would exist. Smaller organisms would be vaguely spherical collections of cells, while larger organisms would be more flattened (if living on land, somehow) blobs of larger collections of cells. This is simply what would result from undirected growth and physics and doesn't interest me. Biological patterns, on the other hand, are complicated and interest me a great deal.

Turing 1952, Fig 3.

Our modern understanding of biological patterns started with a mathematical (and theoretical biology) paper by Alan Turing in 1952. In "The Chemical Basis of Morphogenesis", Turing's key observation is that a system of two components with specific interactions will form semi-periodic waves in space. An "activator" encourages the production of itself and an "inhibitor", which... inhibits the production of the activator. A critical aspect of this reaction is the relative diffusion rates of the components. The activator diffuses poorly, while the inhibitor diffuses readily. All of these characteristics can be described by a small set of differential equations.

Graphic edited from paper.

These differential equations indicate there will be no changes if the space being calculated over is entirely uniform, but the slightest irregularity results in dramatic changes. Regions with slightly too much activator grow and induce inhibition in neighboring regions. Regions with slightly too much inhibitor grow and induce activation in neighboring regions. Quickly the system stabilizes into a series of maxima and minima, visually forming spots/stripes of activation (or inhibition).

Young 1984, Fig 1.

Calculating the change of the activator and inhibitor concentrations over time become computationally expensive if we want to simulate a high-resolution system. Young described a simplification of the activation/inhibition system in his 1984 paper. Because the inhibitor diffuses more readily than the activator, we can think of the neighborhood around each active cell to contain two circular regions. The inner region is where the effect of the activator is dominant. The outer region is where the inhibitor is dominant. These regions can be simplified into two distinct circular regions with uniform activation and inhibition, respectively, instead of the continuous curves described by the differential equations. (Young 1984, Fig 1.)

Young 1984, Fig 2 (left) & Fig 3 (right).

This simplified system is much easier to calculate over a large grid, as a cellular automata. For each pixel, any active cells in the inner region are activating neighbors, while those in the outer region are inhibiting neighbors. If a cell has more activating than inhibiting neighbors, the cell becomes active in the next turn. If a cell has more inhibiting than activating neighbors, the cell becomes inactive in the next turn. If a cell has an equal number of activating and inhibiting neighbors, the cell remains in its current state in the next turn. Depending on how you scale the active/inactive counts, you get a range of patterns. (Young 1984, Fig 2.) Applying the same simplification to an activator/inhibitor model with anisotropic (not the same in each direction) diffusion results in stripes. (Young 1984, Fig 3.)

---

I came across Young's model sometime while I was in high-school, roughly a decade after his paper was published. I was already playing with more conventional cellular automata using software I had written, so I added a module to play around with Young's model.

Fig d1.

One of the first images I produced (Fig d1) illustrates what happens when you bias the ratio of counted active/inactive neighbors. The x-axis biases the active neighbor count (from 0.0 to 1.0), going right. The y-axis biases the inactive neighbor count (from 0.0 to 1.0), going down.

Fig d2.

Next I tried scaling the neighborhood diameters. What I got was seemed to be a random association of diameters to pattern type produced. I didn't clue in to what might be going on until I calculated the effect of smoothly scaled diameters across a single image (Fig d2, top). I realized that the regions I had been using were aliased. Small increases in diameter were resulting in 'random' increases in the number of cells being counted in each region. After rewriting a bunch of code to use pre-defined neighborhood masks, I was able to confirm that using anti-aliased regions produced the pattern I had expected at the start (Fig d2, bottom).

Fig d3.

With the better way to define neighborhood regions, I started adjusting the radius of each region separately. I used some Fourier transforms of the resulting images to help me sort out what was going on, but I didn't save those figures. It turns out that the average feature size in the produced images is the same as the average of the two diameters used in defining the neighborhoods. If the two diameters are further apart, the resulting images show a greater range of feature sizes. By setting the region diameters appropriately, you can design the resulting image to have feature size and feature size variation characteristics of whatever you choose.

Fig d4.

It was also easy at this point to examine the results of anisotropic diffusion in the model. If both regions are squashed in the same direction, the resulting pattern looks more or less normal... only squashed (Fig d4, left). If the regions are squashed in different directions, stripes result (Fig d4, right.). I don't know what happens when anisotropic diffusion is applied to spots vs. stripes. I should explore that.

Fig d5.

Smoothly changing the size of the ellipse neighborhoods produces an appealing network of stripes of different thickness. (Fig d5.) This pattern reminds me of the pattern of stripes [sometimes] seen on the back of cuttlefish.

Fig d6.

Smoothly rotating the ellipses did not produce the results I expected. (Fig d5.) I played around with my code for a while, but never got diagonal stripes to form. I expect the problem is some artifact/bug in the code, rather than something intrinsic about the system, but I still don't know where. It might take rewriting the code from scratch to find my way around it.

Fig d7.

Biological patterns are rarely ever so simple as the outputs of this model have been, so I attempted to simulate added complexity by calculating combinations of Young's model with different factors driving the interaction between one activator/inhibitor pair and another. (Fig d7.) In my example we see two different kinds of spots forming, with interesting interactions between them in the lower middle section.

One of the random interactions I simulated resulted in animated, mobile spots that would wander around the screen and occasionally replicate. Apparently, it doesn't take a very complicated system to start gaining some life-like features. Unfortunately, I hadn't written the program to display the random settings it had chosen. There was no way for me to save them. (I've since then rewritten the code to save any random settings. I will find those amoeba again!)

---

When I started graduate school, I had to set aside my explorations into theoretical biology like this. I'm now done with graduate school and find myself with sufficient free time to post to this blog, so I know I'll find the time to continue exploring the simulation of biological patterns.

Given the wide range of patterns seen in biology, I strongly suspect there are a great many more interesting results to be found exploring combinations of the simplified activator/inhibitor system described by Young's model.


References:

• Turing, A. M. (1952). The Chemical Basis of Morphogenesis. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, Vol 237, No. 640. p37-72.
• www.dna.caltech.edu/courses/cs191/paperscs191/turing.pdf
• Young, D. A. (1984). A local activator-inhibitor model of vertebrate skin patterns. Mathematical Biosciences,72, 51–58.
• www.sciencedirect.com/science/article/pii/0025556484900609
• Wired article
• www.wired.com/2011/02/turing-patterns/
• Graphic source:
• www.sciencemag.org/content/329/5999/1616.figures-only
• My images in a Flickr album
• www.flickr.com/photos/darrenabbey/albums/72157621241644344
• Cellular automata
• en.wikipedia.org/wiki/Cellular_automaton
• Cuttlefish
• www.pbs.org/wgbh/nova/camo/anat-nf.html

Evolution

While thinking about the evolvability of different artificial life simulations, as discussed some in my last posting, I realized that it would be helpful to talk about what is required for a system to evolve. It comes down to four basic traits.

1. Reproduction: Some unit in the system has to reproduce. This unit could be bacterial cells in your gut, or it could be numerical representations in a computer. (Even fire can be described as reproducing when it spreads through a house or forest.)

2. Inheritance: During reproduction, each new unit in the system has to gain traits from its parent(s). The traits could be hidden, as in recessive alleles, or it could be obvious, as in dominant alleles. The number of parents can be one or more than one. (We have two, but maybe some aliens have three or more.)

3. Mutation: At some point in the reproductive cycle, there has to be the potential for changes in the traits (mutations) that are inherited.

4. Death: Death is generally required to remove individuals from a population, thus freeing up room for the next generation. However, there are scenarios where death isn't required. If the population is continuously expanding into new territory, the front-line sub-population can evolve over time without individual death. In this case, the older organisms being left behind fills the same role of actual death.

---

It is relatively easy to prove mathematically that a system with these four traits will experience evolution.

Lets give it a go in a simulation that has a maximum population of four organisms represented by letters and driven by the following rules.

1. Reproduction with inheritance: A -> AA; B -> BB
• A or B can duplicate.
2. Mutation: A -> B.
• A can mutate into B.
3. Death: A -> A 
• Only A can die.

We start the simulation with "A" .

"A" -> "AA" -> "AAAA" -> "AAAB" -> "AAAB" -> "AB" -> "AABB" -> "AABB" -> "ABB" -> "ABBB" -> "ABBB" -> "BBBB"

This may not look like the sort of math you are familiar with, but it is math nonetheless. Math is the manipulation of abstract symbols that represent precise concepts with the extremely rigid rules of logic. 2+2 always equals 4. A system with the described traits will always experience evolution.

Now, this little toy system I've described has an extremely low evolvability. The starting state of the system ("A") does meet the four requirements and thus evolves. However, once the system has reached the final state ("BBBB"), it no longer meets the four requirements and thus cannot evolve further.

---

If you argue that life doesn't evolve, then you are logically arguing that life does not meet one of the four requirements discussed above. Unequivocally, life meets the four requirements.

Life evolves. The math doesn't provide any other possibility.