Virtual Chemistry

If you’re interested, Virtual Chemistry can now be played with or downloaded here.

A little while ago, I was, once again, idly messing around in NetLogo, having fun with the layout-spring method, which basically treats agents like the nodes of a network of springs, and the links between them as the springs themselves, then computes the dynamics and changes the positions of hte nodes accordingly. Eventually, inspiration struck, and I wondered if it would be possible to build a simulation of chemical bonding. I spent about an hour writing and revising the code, and thus, Virtual Chemistry was born.

It works like this: at the beginning of each simulation run, a random number is generated. This is the number of unique elements that will appear. Arrays are created that store the properties of these elements, namely: complements (which tells the element which other elements it can form bonds with), and maxbonds (which, obviously, tells the element the maximum number of bonds that it can form). Every step, an atom checks for other atoms within the radius “interact-radius” (which the user can set using a slider). If any of the nearby atoms are in the atom’s complement list, there is a ten percent chance that it will form a bond, as long as it hasn’t reached the maximum number of bonds. Atoms also wiggle around randomly to simulate temperature effects.

It took me a while to get this to work and to get rid of a few irritating bugs, but eventually, I got it working, and, not to be immodest, but I was impressed with the kinds of behavior it could produce.

As I watched the simulations unfold, I observed behavior that I hadn’t even considered when I was building the simulation. For example, when the simulated space was packed to a sufficient density with atoms, it began to experience real pressure effects. The way the physics works, bonds hold atoms together, but atoms repel each other. At the start, with the atoms at a low density (and a low enough temperature), certain “molecules” formed. As I added atoms (effectively increasing the pressure), some structures that had been unstable became stable, and some of the structures that had been stable destabilized. I was fairly happy, because I realized that what I was seeing were actual phase changes in my virtual material. Similar phase changes occurred when I changed the temperature, too.

One experimental run in particular is very illustrative. I call it Triad World. In this world, there is only one element (“Element 0”), and each atom of Element 0 can have at most three bonds. Even in a simple world like this, I observed some interesting things.

Here, you can see the atoms, and the three-atom triangles that give Triad World its name. As you can see, these dominate, but there are also other “molecules” which occur quite frequently. The scene above arose at a low temperature and a fairly low pressure.

After adding a few atoms (and thus increasing the pressure), new stable structures emerged. Note the large chain near the center of the image. Only the restraining force exerted by the surrounding atoms prevents the heat from tearing this chain apart. The same force, though, also prevents the two-atom dimers from forming triad,s which they normally would have done.

There are a lot of potential applications for a program like this (more so, perhaps, than my zombie infection simulator), and it’s certainly a lot of fun to play with. The problem is, running the simulation with more than about two hundred atoms slows it down pretty badly. But in order to simulate anything as complicated as, say, a rudimentary biomolecule, thousands of atoms or more would be required. With that in mind, here’s my to-do list of improvements:

  • Optimize the code, if possible, to make it run faster.
  • Write a better bond-forming routine. I’m not happy with the one that’s currently implemented.
  • Tune the attraction and repulsion strengths. Right now, I get the vague feeling that the weird interatomic forces that exist in the simulation are preventing some interesting structures from forming. I know for a fact that they’re preventing any but the most rudimentary chemical reactions.
  • Implement some kind of energy system. Right now, the only energy comes in the form of the random motions induced by the heat. I’d like to make it so that forming bonds consumes energy and breaking bonds releases it.

That’s all for now. I’m hoping to have the complete program up on the NetLogo website soon, and when it’s up, I’ll publish the link.

SimHeart 3.0

Some months ago, I wrote a series of posts about the behavior of my NetLogo heart simulation (namely, this one and this other one). SimHeart has been on the backburner since then, partly because I was making some effort (some) to pay attention to my classes, but mostly because I didn’t have any new ideas to implement. Then, a couple of weeks ago, I was messing around in NetLogo, trying to figure out how to get the cells of the grid to obey a switching-on rule similar to the kind of rules you find in neural networks. That, as it turned out, was a huge bust, but I re-used some of the code to build SimHeart 3.0.

This is a pretty radical revision, but, behaviorally, it’s pretty much the same as the old simulator. Still, there are some changes:

  • The cells, not the agents, do most of the work. Whether they fire or not is determined by a stochastic (random) sigmoid function. If a lot of a particular cell’s neighbors have high potential (the primary variable; actually, three varaibles, one for each system, but I’ll get to that later) and the cell’s potential is low enough, then the cell has a very high probability of switching its potential to 1. If not, the potential steadily drops.
  • The AV node is now fully simulated. I mentioned in my previous posts that I wanted to do this, and since the atria and the ventricles are simulated as two different kinds of potential, all I had to do to simulate the AV node was to introduce a third variant of the potential variable. Now, a heartbeat starts at the SA node (meaning a particular cell’s potential is set to 1), the wave travels until it reaches the entry cell for the AV node, triggers another wave (one that’s morphologically different, which I’ll talk about in the next bullet point) in the AV node, which travels to the exit cell, triggering the ventricular beat. This solution has the nice feature of taking care of the AV node’s natural delay for me, as well as realistically limiting the heart rate that the AV node can transmit from the atria to the ventricles. This is, by far, the most important addition.
  • The model now runs faster. Not much faster, but every bit counts, and given how much simpler the code is now, the model is now much more compact and elegant.
  • The model is more realistic. The stochastic potential-based cells are a lot more like how a real heart works than was the old model.
  • The model now takes into account the different sizes of the heart’s components. Because of the way the model works, the longer a cell’s potential takes to decay, the smaller of the particular component being simulated. Thus, the potential in the AV node takes a very long time to decay, making nodal arrhythmias essentially impossible (and making the AV node behave as though it’s very small); the potential in the atria takes a moderate amount of time to decay, making atrial arrhythmias much less common (and making the atria behave as though they’re larger than the AV node, but smaller than the ventricles. I think you see the pattern here.), and the potential in the ventricles decays quickly, making them fairly arrhythmia-prone.
  • The arrhythmias have changed. This isn’t necessarily a good thing, but it’s not crippling, either. The model now produces ventricular tachycardia much more readily than it did before, but the incidence of ventricular fibrillation has been reduced proportionately, and can be hard to differentiate between tachycardia. Unless the parameters are adjusted, atrial arrhythmais don’t occur at all. Still, the model can now incorporate junctional tachycardia (where a spiralling, self-sustaining wave in the AV node stimulates the ventricles to beat too quickly) can now be simulated easily.

And, not mentioned on the list above, the model is now prettier, too, since I was apparently experiencing some aesthetic inspiration when I was writing it. So, in SimHeart tradition (can I call it a tradition after three posts?), I present: the screenshots (note: ignore the weird screwed-up areas. They’re just explanatory text that apparently didn’t transfer into the screenshot. Don’t worry, you’re not missing anything, it’s all explained in the text below):

Periodically (with the period determined by the “s-rate” slider next to the view), the cell labeled “SA” sets its potential to 1, generating the atrial pulse (the red wave). If, for some reason (for example, if the atria are in fibrillation and there’s already a wave passing through the SA node), no wave will be generated.

When the pulse passing through the atria hits the cell labeled “AV1,” if there’s not a wave passing through it at the moment, it triggers a pulse in the AV node (the green wave).

If all goes well, the AV pulse triggers the cell labeled “AV2” which activates the ventricular pulse (the blue wave).

Of course, all does not always go well:

In this case, something has gone wrong with the electrical wave as it moves through the ventricles, causing a deadly spiral wave to form. This is simulated ventricular tachycardia, which will (probably) eventually decay into ventricular fibrillation. As you can see on the ECG, the atrial pulses are still trying to get through, but the AV node can’t activate the ventricles. Although it’s not pictured here, a quick press of the “Defibrillate” button took care of the arrhythmia.

Since the model’s now so much simpler, it’s easier to simulate diseases and treatments. Note the row of blue buttons on the right side of the images. There is, of course, the defibrillate button, which is remarkably effective at terminating arrhythmias. Below that are the “Increase Adrenalin” and “Increase Antiarrhythmics” buttons, the first of which makes the heart depolarize faster, and the second of which makes it depolarize more slowly. The first, as you might expect, makes the heart more arrhythmia-prone, and the second, obviously, makes it less arrhythmia-prone.

The set of buttons below that simulate various cardiac illnesses. “AV Block” makes the AV node’s cells depolarize very slowly, meaning that, when the sinus rate is set to a very high value, not all of the beats are able to get through. This, if I’m not mistaken, represents 2nd-degree AV block, type 1. The button below that, “Sinus Tachycardia”, sets the sinus rate to a very high value, simulating either a disease process or the effect of very strenuous exercise or other stress on the heart. “Long QT Syndrome” is a sort of bastardized version of the real disease, but has the same effect, making the ventricles dangerously more prone to arrhythmia.

All in all, I’m far happier with the new version than I was with either of the old versions. There are still some improvements to be made, however, and I’ll post updates as needed. Soon, I hope to have the model up on the NetLogo website, so that everybody can fiddle around with it.