Category Archives: MatLab
“I am a spec on a spec of a planet in the spec of a solar system in one galaxy of millions in the known Universe. My Universe is like a ring in the desert compared to the Footstool of God. And the Footstool like a ring in the desert compared to the Throne of God.”
The amount of research done on brains, neurons, neuro-chemicals, etc….is astounding. No matter whom your favorite neuroscientist is you can rest assured that his or her knowledge barely scratches the surface of all that is known about human brains, or even Bumble Bee brains for that matter. But even all that’s known is but a “ring in the desert” compared to all that can be known–to all that will one day be known.
Hodgkin and Huxley developed their model of axonal behavior in the European Giant Squid in the late 1940s and published it in 1952. It was a game changer for two main reasons–first because it accounted for the action of Sodium and Potassium channels, and second because it stood up to experiment. It is a shining example of mathematical biological modeling from first principles.
FitzHugh and Nagumo based their model partially on the Hodgkin-Huxley equation, partly on the van der Pol equation, and partly on experimental data. While it’s not as accurate as the H-H equation relative to experiment, it’s much easier to analyze mathematically and captures most of the important attributes. It’s also been successfully applied to the behavior of neural ganglia and cardiac muscle. There are many different interpretations of the FitzHugh Nagumo systems of equations. The readiest on-line tool for examining them was created by Mike Martin.
In this post I’m working with a slightly modified version of the forms described by Ermentrout and Terman:
As you can see, this is a system of ordinary differential equations. The second equation has the quality of stabilizing the first. In a sense, it acts like an independent damper on changes in voltage when current is applied to an axon. Like the van der Pol equations, it’s non-linear and deterministic, and doesn’t easily lend itself to analytical solutions.
One thing you can’t see with Dr Martin’s tool is an interesting phenomenon that frequently occurs with these sort of equations: the Hopf Bifurcation. The only equation parameter that can be easily changed during experiment is the applied current. By programming the equations and calculating across a range of currents, it can be determined when an applied current produces unstable voltage oscillations. The point at which it becomes unstable is known as the critical point or Hopf point, and as long as the region of instability doesn’t extend to infinity, there will be one on each side. According to Steven Baer, the critical points for Ermentrout and Terman’s system are found at:
Enough about all that. Using parameter values of a = 0.8, epsilon = 0.5, and gamma = 0.2, I punched it out using my MatLab program FHNBifurc, and also ran it with XPP. XPP is no-frills and a little bit glitchy, but it’s awesome for solving systems of ODEs and it’s a champ of a program for bifurcation analysis. The installation is a little tricky, but it’s totally free–no spyware or anything else attached. If you haven’t at least toyed around with it, you should.
The MatLab program works by drawing the diagram in two directions–from the right in red and from the left in blue. The calculated critical points are shown as white stars. Ideally, the bifurcations should begin at those points, but as you can see, they don’t show up perfectly. That’s the image on the left. The best use of the program I’ve found is that it can be used to search across a wide range of currents to determine if and where instability occurs. I ran the same stuff with XPP’/Auto, and it gave me the figure on the right. Auto’s kind of a bear to get to the first time, but do it once and you’re set for life.
I also did a final run with Matlab, this time a smaller range with a lot more iterations. Took my computer about 10 minutes to complete it, so make sure you’ve got enough memory and processing speed before you try it. I edited-in the critical point solutions, FYI:
You can see why the FitzHugh-Nagumo equations are called eliptical. The MatLab program is great for scanning a wide range of values and locating the range of the bifurcation. Auto is way primo for drawing their shapes. Here’s the code for the XPP file:
# FitzHugh-Nagumo equations # fhnbifurc.ode dv/dt=-v*(v-1)*(v-a)-w+i dw/dt=e*(v-g*w) parameter a=0.8, e=0.5, g=0.2, i=0 @ total=50000, dt=1, xhi=50000., MAXSTOR=100000,meth=gear, tol=.01 done
The MatLab call looks like this:
a = value of a eps = value of epsilon gamma = value of gamma
I1 = initial current value IF = final current value
tspan = time span over which current is scanned
tstep = how short the time steps should be for the calculation
The specific call for the 0 to 10 scan was:
and for the more detailed diagram:
As always, the complete code is below. Cheers!
PS–Did you like The Sexy Universe on Facebook yet?
function FHNBifurc(a,eps,gamma,I1,IF,tspan,tstep) steps = tspan/tstep; Iramp = (IF-I1)/steps; V(1) = 0; W(1) = 0; I(1) = I1; for n = 1:(steps-1) In = I(n); Vn = V(n); Wn = W(n); fV = Vn*(Vn-1)*(Vn-a); dVdt = -fV - Wn + In; dWdt = eps*(Vn - gamma*Wn); V(n+1) = Vn + dVdt*tstep; W(n+1) = Wn + dWdt*tstep; I(n+1) = In + Iramp; end IFwd = I; VFwd = V; WFwd = W; V(1) = 0; W(1) = 0; I(1) = IF; Iramp = -Iramp; for n = 1:(steps-1) In = I(n); Vn = V(n); Wn = W(n); fV = Vn*(Vn-1)*(Vn-a); dVdt = -fV - Wn + In; dWdt = eps*(Vn - gamma*Wn); V(n+1) = Vn + dVdt*tstep; W(n+1) = Wn + dWdt*tstep; I(n+1) = In + Iramp; end IRev = I; VRev = V; WRev = W; Vcr(1) = (a+1-sqrt(a^2 - a + (1-3*eps*gamma)))/3; Vcr(2) = (a+1+sqrt(a^2 - a + (1-3*eps*gamma)))/3; for n=1:2 Icr(n) = (Vcr(n)/gamma) + Vcr(n)*(Vcr(n)-1)*(Vcr(n)-a); end Vcr Icr plot(IFwd,VFwd,'-r',IRev,VRev,'-b',Icr,Vcr,'*w') xlabel('Current (I)') ylabel('Voltage (V)')
This was originally posted at The Cameron Hoppe Project. It’s popularity inspired this site. It’s been updated since.
“And He will raise you up on eagle’s wings,
bear you on the breath of dawn,
make you to shine like the sun,
and hold you in the palm of His hand.”
Felix Baumgartner completed his much pre-hyped jump from 128,100 feet, achieving a top speed 1.24 times the speed of sound. It must have been an amazing ride, with the bill footed by Red Bull and everything. I am neon green with envy. Sure, it could be death to try, but it would be worth it just to see the blue-glowing world one time. Besides, the life insurance is paid.
Gravity creates constant downward acceleration. Friction with air produces drag that runs the opposite direction in which one is moving, and is proportional to velocity squared. So the sum of all the forces on Felix, as with any skydiver, during his fall were equal to his mass times his actual acceleration. Since drag is proportional to the square of velocity, the drag force and the gravitational pull become equal, and acceleration reaches zero. This is known as terminal velocity. Not nearly as exciting as the name sounds. I know; I was disappointed, too. In equation form it looks like this:
The drag coefficient includes half the surface area of the guy in the suit and a proportionality constant. Based on his reported time in free fall and an area of 4.3354 square meters I found this to be 1.15, which is actually quite reasonable.
The next issue to deal with in the problem is the air density. Anybody who’s ever been up a tall mountain or a plane ride knows air gets thinner with elevation. As a result, drag forces are higher near the ground than at the elevation Baumgartner jumped from. Without that included, we’re only left with the beautiful sky, without any interesting math. So we can expect that Felix’s speed went way up, and then it actually decreased until he pulled his parachute at 270 seconds. Then it really dropped. The air density equation is described below:
There’s a ton more I just didn’t have time to do. For example the gravitational force also decreases as elevation increases and the temperature gradient is probably not a constant from the ground to the stratosphere, etc, etc. But time is short. Well, we’ve got two coupled differential equations. There’s only one thing to do–code it and find out how it looks! I punched this one out in Matlab:
The scale of this is in tens of thousands of feet, so you can see where he pulls the cord at 270 seconds and about 8400 feet. At that point, his descent slows way down. Now, the velocity of the fall:
Two important spots here. On the right, you can see where his chute deploys. On the left, you can see where he reached terminal velocity in the upper atmosphere about 50 seconds in. This was his maximum velocity. After that, the rising density of the atmosphere continually slowed him down. Kinda like being married.
You can see the model is a little off; the maximum velocity should be 1223.75 feet per second, while I’ve got him maxing out around 1175. It’s probably due to my crude modeling of gravitational and temperature gradients at high elevation. What can I say? There’s only so many hours in a Sunday afternoon. The model is named for Baumgartner but it can be used generically for any object falling through a gas. It works out everything in SI units, then converts to feet at the end before plotting. Changing the plot command to small letters puts it in SI. The rest is pretty straightforward. The function call is:
r = Radius of object L = Length of object CD = Drag Coefficient
p0 = initial position u0 = initial velocity a0 = initial acceleration
rhoR = density of object rho0 = density of air at sea level
T0 = Temperature at sea level in Celsius Tf = Temperature at p0
The call I used was:
PS–Are you following @CameronHoppe on Twitter yet? If you follow me, I’ll follow you. Plus you should leave a comment. Just sayin’.
function Baumgart(r,L,CD,p0,u0,a0,rhoR,rho0,T0,Tf) tstep = 0.00001; Tgrad = (Tf - T0)/30000000; for n=1:30000001 t(n) = (n-1)*tstep; u(n) = 0; p(n) = 0; end n = 0; for n=1:27000001 a(n) = 0; rho(n) = 0; T(n) = 273.15 + T0 + (n-1)*(Tgrad); end a(1) = a0; u(1) = u0; p(1) = p0; rho(1) = rho0*exp(-0.0284877*9.8*p0/(T(1)*8.3144621)); n = 0; V = pi*L*r^2; A = 2*pi*r*(r+L); g = -9.8; k = 0.5*CD*A; massR = rhoR*V; for n=2:27000001 u(n) = u(n-1) + tstep*a(n-1); p(n) = p(n-1) + tstep*u(n); rho(n) = rho0*exp(-0.0284877*9.8*p(n)/(T(n)*8.3144621)); a(n) = (massR^-1)*(V*(rhoR-rho(n))*g + k*rho(n)*(u(n)^2)); end uS = u(27000001)*.25; for n = 27000002:30000001 u(n) = uS; p(n) = p(n-1) + tstep*u(n); end for n = 1:30000001 P(n) = p(n)*3.28084; U(n) = u(n)*3.28084; end A k m = plot(t,P); xlabel('Time (seconds)') ylabel('Position (feet)') set(m,'LineWidth',2)