Concerning Contagion

A measles outbreak has been in the news lately, and since the disease was previously considered eradicated in the US, I wanted to try a simple model of disease spread with and without vaccination. To represent a population, I used Python's NetworkX package, which gives us nodes (people) connected by edges (contact that can spread disease). In addition to vaccination status, I wanted to see how different levels of contact affected spread, so I generated the layouts using the Watts-Strogatz Graph, which takes a given number of nodes n, a community size k, and a rewiring parameter pwire. The nodes are arranged in a ring, with each connected to its k nearest neighbors. Then we randomly switch some of those connections according to pwire to randomize the interactions a bit.

Turning now to the nodes of these networks, the model for people I chose tracked a few qualities: Current infection, vaccination/immunity status, susceptibility to infection, and infectious period. Based on some stats I found from the CDC, I decided after being infected a person would be sick for 22 days, but only contagious for days 10-18. While contagious, each day a person's neighbors can become infected according to their susceptibility (0.5% for immune, 50% otherwise). I started the model with 10% of people infected, and ran until all cases were gone or a maximum of 300 steps (days). I ran several different values for k and pwire, as well as different vaccination levels.

After running the simulations, I tried to come up with some summary statistics. The first thing that came to mind was the fraction of people that get sick:

The x-axis shows the k parameter for each run, which represents the number of people each person has contact with. For the lower values of k, this shows clear separation in the different levels of vaccination – Low vaccination leads to high infection rates. With enough connection though, everyone encounters an unvaccinated person and gets sick.

That's just the number of people who ever get sick, but we can also look at how long the disease spreads before being eradicated:

Again we see a big uptick for high community size, to the point where disease is still spreading after 300 steps. For the smallest size, we again see separation in the vaccination rates, but there's also an effect from the pwire parameter. This is easier to understand if we look at a map of one of these communities:

Blue dots are people who have never been ill, red are currently sick, and green are vaccinated, or recovered. The twisting shape shows how the pwire parameter has connected a few people on opposite sides of the initial ring. The black dot in the plot above shows that for low pwire, the disease has to make its way around the whole ring, increasing the time to eradication. We can also look at a case with larger communities and more interconnection:

In this case, almost everyone immediately gets sick, which you might imagine would reduce the time before the disease is eradicated, but the increased connection also results in many more breakthrough infections.

I'm really impressed with the complex behavior this simple model was able to show! That said, I fully recognize it's just a toy, and should not be used for setting policy. I'm not headed to that farm upstate just yet...

Make Like a Tree and Get Out of Here

Swarthmore's Scott Arboretum traditionally gives incoming students a plant to care for in their dorms. Shockingly, the Hawaiian Schefflera I received in 2007 is still going 17 years later! Over the years as I've moved from place to place and kept the plant in different environments, I've been impressed by its ability to track the sunlight, frequently growing lopsided as it reaches toward the nearest window until I think to turn it. This has resulted in some twisting, gnarled branches:

If you've been reading this blog, you can probably guess where this is going – I was curious if I could make a simulation of my plant's heliotropic tendencies. I decided to model the plant as a collection of connected branch segments, each a fixed length and pointing at an angle relative to the vertical. At each step, we iterate over all the segments and pick an action:

• If the segment has no children, i.e. it's at the tip of a branch, we add a new segment on the end with a probability p_grow/size, where size is the number of existing segments.
• If the segment does have children, we add a new one with probability p_sprout/size.
• If neither of those occur, we adjust the angle of the branch to point closer to the sun's current position. The adjustment is proportional to how far off the angle is, how many branches are on the end of this one, and a constant stiffness for the plant.

I tried a bunch of values for the different parameters until I landed on a range that gave plants looking reasonably similar to the real thing (click to enlarge):

The numbers along the top give the stiffnesses, and the ones on the left give the sprout probability. The sun moves back and forth sinusoidally, which you can see in the snaking of the plants. I wasn't able to get my digital plants to spread as much as the analog one, possibly because I'm not accounting for the plant casting shadow on itself, but I'm still pleased with the results – The top center one seems particularly good. If you'd like to try for yourself, the code is here.

Stopping Traffic

Since coming back to Florida and commuting every day, we've noticed that the traffic lights stay on a particular direction for a long time before switching. This is especially annoying as a pedestrian (who dislikes jaywalking), since it means the walk from the parking lot to my office can vary a lot depending on what lights I hit. I was curious whether I could design a simulation to test the effect of the light duration on the time cars spend waiting.

The setup is a 4-way intersection with lights that swap red/green with separate green durations for North/South and East/West. Each road has a rate cars arrive, and each car randomly chooses to go left, right, or straight. As cars pass through the intersection, we can keep track of how long they had to wait since arriving. Below is an example simulation – The arrows are the lead cars on each road, pointing in their planned direction and colored to indicate how long they've been waiting.

We've got a lot of variables going on here, so it took me a while to figure out how to assess the results. It would take too long to try every combination, so I just sampled from a reasonable-sounding range. I considered using a corner plot to see the effects of the inputs, but the sampling was too sparse to see a trend. One way to reduce the number of parameters is to consider where each car is coming from, and use the traffic rate and light duration from that side. Then I thought about the units: the traffic rate measures cars/time, and the light duration is a time, so multiplying them gives an average number of cars during a green light. Plotting this with the waiting time shows a clear trend:

I separated out the different turning directions, expecting that left turns would need to wait longer, but that doesn't seem to be the case. This shows that to minimize the waiting time, we need to keep the rate*duration small, i.e. for busy roads the lights should change rapidly, which is not what I expected.

This also implies that the Florida lights are expecting a low rate of traffic – Not our experience the past few weeks coming home during rush hour and hitting gridlock! Of course, this model may be inaccurate, and I should stick to gravitational waves, rather than civil engineering.

Head's Up

Recently, I saw my father-in-law Scott pour a bottle of beer into a glass, and I was fascinated by the relationship between the rising beer at the bottom, and the foam moving on top. I was curious if I could model the dynamics involved, so I decided to check if anyone else has tackled the problem. I found an article from a brewer discussing some of the steps involved, and I decided to split the process into 4 parts:

1. For each bubble in the foam, apply forces from neighboring bubbles, and gravity pulling down.
2. Drain liquid from higher to lower bubbles, based on the content of each. Bubbles at the bottom drain into the liquid beer.
3. If neighboring bubbles each have low liquid content, merge them into a single larger bubble.
4. Bubbles with low liquid and/or large size pop, adding their liquid to the beer at the bottom.

To represent these bubbles, we need to use sphere-packing to figure out how they fit in the glass. I've talked about the concept before, but since we're making a simulation this time, I found the package spack for Python. This handles keeping track of where the bubbles are, and can draw them for us. It also calculates the forces between the bubbles, based on their separation and radii.

To apply these forces, we need the mass of the bubbles – They're made of gas and liquid, but the liquid mass will far outweigh the gas. We're already keeping track of liquid content for the other steps, so we can use that for the mass and then displace each bubble based on the total force divided by its mass.

The spack package keeps track of which bubbles are adjacent – For each pair, we can drain liquid from the upper to the lower. The proportion I settled on was that at each step the lower bubble would get 51% of the total moisture, and the top one would be left with 49%. For bubbles with nothing below them, they drain the full amount to the liquid at the bottom. Similarly, we can use the adjacent list to pick bubbles to merge – Based on the article I linked above, I decided their circumferences would add, rather than areas. We select which bubbles merge based on their liquid content – Dryer bubbles merge more easily. Big, dry bubbles can also pop, giving up their liquid to the beer at the bottom.

We start off the simulation by filling the glass with a bunch of small bubbles, then let the rules outlined take over. I ended up using 1000 bubbles to start with, which takes a bit of time to get through, but I'm impressed with the results:

We can also measure some averages over the course of the simulation. First, we can simply count the number of bubbles:

The rate is fairly constant, despite the various dynamics going into determining the merging/popping. We can also look at the average size of bubbles:

I was a bit surprised by the sudden rise in the average size at the end, but I think it may be related to the merging of the bubbles: They get larger and fewer, resulting in a compounded effect on the average, and near the end, we have far fewer bubbles, making the average more sensitive to change. The average liquid content has a similar knee at the later times:

I have no idea how accurate this model is, but it does seem to follow the events outlined in the article: Bubbles merge, dry out, and pop, resulting in a shrinking head of foam, and growing reservoir of liquid at the bottom. Clearly I need to gather more data – Cheers!

Dirty, Disgusting, Filthy, Lice-Ridden Boids

[Title from The Producers]

Long ago I had a screen saver with a simulation of bird flocking behavior – A group of 2000s-era 3D blocks would fly in formation, land, and take off. I was recently reminded of it, and grew curious how it was made. My best guess is that it used a model developed in 1987 called Boids. The model consists of a group of agents (boids) that each act according to a set of simple rules. In this case, the rules are

• Separation: Avoid flying into nearby flockmates
• Alignment: Fly in the same average direction as nearby flockmates
• Cohesion: Fly toward the average position of nearby flockmates

Along with these, I also included a target location that all boids want to get to. There's a lot of details unspecified here, like what "nearby" means, and how these various rules are weighted, e.g. cohesion and separation can be directly opposed at times. I looked at a couple implementations I found online, and tweaked my own model until I got reasonable looking results.

In the simulation below, the boids are all heading for the center of the map, and "nearby" covers about 1/4 of the map. I used "toroidal boundary conditions" which means if a boid goes off one edge, it wraps around to the opposite side. You can see that a few boids sometimes break out of the "nearby" region and head off on their own before rejoining:

I was pretty happy with how this worked out, but I felt like I should give it more of a physics spin. It occurred to me: What if these were relativistic boids? If they're moving at a significant fraction of the speed of light, then the observations they use to follow the rules above are based on the light arrival time, which may not reflect the current motion of a particular boid. What this amounts to is that based on the distance between two boids, their observations are delayed by a certain time. Here's what happens in that situation:

I was really surprised by the results – I expected that they would have a much harder time sticking together, since their actions would be too imprecise, but instead they show better clustering. I think this is because by staying closer, they get more up-to-date measurements, so it ends up being a positive-reinforcement.

I figured I could get a more qualitative comparison by computing some summary statistics. First I looked at the average distance between the boids at each time:

This shows the relativistic boids bunching up quickly, but then oscillating sharply around a higher value than then lowest the non-relativistic boids are able to hit. We can also look at the deviations in the boid's headings for the two cases:

This shows the relativistic boids frequently going in different directions, which we can see in the animation above: The tight grouping requires quickly changing direction to stay in one place.

I'm not sure whether any of the space organizations have plans for high-speed probe swarms, but if they do, I hope this will serve as my grant proposal!

Self-Jamming Cars

This week I wanted to try another Complex Systems-style simulation, this time based on something I had originally seen on Mythbusters, but others have tried with similar results: Traffic jams that appear out of nowhere simply due to drivers varying their speed. The system is fairly simple: Some number of cars drive on a circular track, trying to go as fast as possible while maintaining a safe stopping distance and obeying the speed limit. However, they may not all accelerate at the same speed, and can brake unexpectedly. These error factors result in some interesting effects, including waves of slow speed that travel backwards around the track.

To simplify things, I assumed the cars were either accelerating at maximum, or braking at maximum. They would decide based on the stopping distance from the car ahead of them:

where v is the car's velocity, and a_b is the braking acceleration. If the distance to the next car is less than this, we brake, otherwise we continue accelerating.

The controls you'll find below are the number of cars on the track, the maximum speed they'll go, the rate of acceleration, the rate of braking, the size of the variation in acceleration rate, and the rate at which that variation changes. This last factor is needed because if we change the acceleration error at every step, it tends to average out and have little effect. There's also a brake button that makes the red car slow while holding it down. As with many of these Complex Systems topics, I'm always surprised by the dynamics that emerge with just a few simple rules. Be sure to post a comment if you find something particularly weird!

Mind Your Pizz’s & Queues

Most Friday nights, Marika and I like ordering pizza for dinner from our favorite Gainesville spot, Satchel's Pizza. Being a Friday night though, many people have the same idea, and we can often get stuck with long wait times. This past Friday, we decided to get our order in earlier, ended up getting in before the rush, and our pizza was ready sooner than we would have liked. This situation, strategizing our moves relative to others trying to do the same thing, reminded me of a project I was part of when I first started at the University of Michigan. It was in a parallel/sub-field to Physics called Complex Systems, where we try to apply techniques from Physics to other systems. In this case, we were looking at co-adaptation and co-evolution. Our model involved simulating a group of agents that would act on an environment with different strategies. Based on the strategy, they would change the environment, and get some benefit or penalty, then try to change their strategy for more benefit.

How does this relate to pizza ordering? The other customers and I are the agents, acting on the pizza place. We choose when to place our orders, and based on how many orders are ahead of us, it takes the pizza place a certain amount of time to prepare the order. We then change our ordering time according to how close to dinnertime our pizza finished. The adjustment I settled on was

This says we adjust our ordering time by 10% of the time between when we want to have dinner, and when we actually got our pizza. I put the 10% in there to make the transitions a bit smoother.

The pizza place I modeled as a set of ovens and a queue. When a customer orders, they're added to the queue, and each empty oven will serve the first customer in the queue, and take 20 minutes to cook the pizza before becoming empty again. I was curious how the customers would change their order times, and what kind of wait times would be involved. I decided to run the simulation with a couple different numbers of ovens, and see how things changed. I decided all the customers would try to eat at time t=0, and use t=-20 minutes as their first order time. Here are the average order times used by the customers at each iteration:

After some initial transients, each of the oven cases settles into a sinusoidal pattern, but with different amplitude and frequency. We can also look at how far off the customer was from their dinnertime:

A couple weeks ago was Homecoming weekend, and the wait times did indeed get into the 2 hour range! Where these results get really interesting is if we combine the previous two plots into one showing the relationship between the order time and the wait time:

This has the appearance of an attractor, a common feature of complex systems, where the state will trend toward an equilibrium point, but not necessarily reach it. This is the average over all the customers though, so we can also look at the individuals over time:

It seems the other Satchel's customers and I are trapped forever in a cycle of never getting our pizza quite when we want it!