1. ** Finding Roots Algorithms **

2.** Ordinary Differential Equations **

3.** Using the Runge-Kutta Method to solve the Hodgkin-Huxley Model **

__The Problem__

I choose to describe a problem discussed in an undergraduate course I have taken in the EE department - Introduction to Biological Signals and Systems.

In physiology, an action potential is a short-lasting event in which the electrical membrane potential of a cell rapidly rises and falls, following a consistent trajectory. The shape of a typical action potential is depicted in the following figure. The membrane potential remains near a baseline level until at some point in time it abruptly spikes upward and then rapidly falls.

The Hodgkinâ€“Huxley model, or conductance-based model, is a mathematical model that describes how action potentials in neurons are initiated and propagated, through a system of four first-order non-linear differential equations:

where V is the membrane potential, I is the externally applied current. V_{Na},V_{K},V_{Cl} are the ions' resting potential. The current of each ions type is proportional to the potential difference between the membrane potential and the rest potential. The parameters g_{i},α_{k},β_{k} and the rest potentials are dependent on the neuron type and are found experimentally. Following is an example set of parameters:

We wish to solve the above set of equations for V(t), given a set of initial conditions and a given externally applied current.

__Algorithm__

I would chose to solve this problem using the fourth order Runge-Kutta method. This method is chosen because it is more accurate than the Euler method and since the equations are not linear, the implicit methods can be difficult. The steps for solving would be the following:

For a given set of initial conditions and not too long a time vector, this should not be too long of a calculation, and I would plot V(t) at the end. But, if we wish to examine the response to changing external currents for a long period of time, I might choose to plot the progress of V(t) at every iteration.

__Accuracy and Validation__

The accuracy is the order of dt^{4}, so the smaller dt the better. But, as we've learned, this can be a problem if dt is smaller than the precision of the arithmetic in the computer.

Another problem might occur if the initial conditions and external current are such that will create chaos. Reduction of the above equations to two dimensions under certain assumptions can be examined in state-space representations and there isn't always a steady state. In real life this cannot happen (since n and h bring V down) but if the given conditions are not realistic one might receive chaos. This should not happen.

__Question 1.5:__I chose to write the program in C. The code can be found here.

I investigated the behaviour of the Secant method with changes in the initial guesses of the root, and the variety of cases that were tested are described in the following table:

It can be noticed, that as long as one of the guesses is close to the root, there is pretty fast convergence, and that the closer the guess to the root, there is a faster convergence. Also, it can be noticed, that the value of x1 is slightly more significant than the value of x0. But more importantly, since this is a simple quadratic function, no matter how large the initial guess is, the algorithm converges to the correct answer. I even tried running it with x_{0}=1000001 and x_{1}=1000000 and received the answer within 32 iterations.__Question 1.6:__

This program is also written in C. The code can be found here.

When I ran the code, I saw that for initial guesses under 1 or very close to 1 the algorithm converges, but even for an initial guess of x=2 it diverges.

For x>>1 the method does not converge. This can be easily understood from the tanh graph - for x>>1 tanh(x) is approximately constant, meaning f'(x)=0 and therefore the expression f/f' goes to infinity.

Now we wish to derive the critical value of the initial guess above which convergence will not occur.

This will happen if x_{i+1}is farther away from zero than x_{i}, meaning: |x_{i+1}|>|x_{i}|. Let's define g(x)=f(x)/f'(x), then we shall demand |x_{i}-g(x)|>|x_{i}|. The solution is g(x)>2x. Now, if x>x_{0}, this directly implies that x_{2}<-x_{1}and x_{1}< -x_{2}< x_{3}etc. The NR algorithm diverges. In our case f(x)=tanh(x) and f'(x)=sec^{2}(x), so the critical value of the initial guess is received by tanh(x)/sec^{2}(x)=2x and the solution is x_{0}=1.08866. And indeed, for x_{0}=1.08866 the code produces an answer of divergence but for x_{0}=1.08865 we get the correct zero in 13 iterations.

I tried solving with the secant method as well (in the same code) and found that it converges for slightly higher initial guesses, but also has a limit.

For an initial guess of zero, the NR method converges in a single iteration. For the secant method, if x1 is zero it converges in a single iteration, not matter how large x0 is, but if x0 is zero and x1 is not it will take more than one iteration, depending on the x1 guess.

2.

__Question 2.8:__

I changed the potential to a central one: V(x,y)=2*sqrt(x^{2}+y^{2})^{5}. Following is the result for E=.1 and initial conditions y_{0}=50, p_{y0}=2. Please excuse me for not changing the titles of the graphs, I didn't want to mess with the code too much.

I tested this for many different initial conditions and found that if y_{0}is smaller the radius of the cross section of the toroid is bigger (infinity fr y_{0}=0). For example, following is the same potential for y_{0}=20 and the same p_{y0}=2.

Note that the radius is larger than before.

But, for the momentum it's the opposite - the larger p_{y0}the larger the radius. For example, following is the same potential for py_{y0}=20 and the same y_{0}=50.

In conclusion, the shape of the trajectory in the phase space does not change, but the initial conditions affect only the geometry a little. We recieve the toroidal shape, just as expected from Fig. 2.2 in the text book.

__Question 2.10:__

For E=0.025 with time stamp 0.1 and initial conditions y_{0}=2, p_{y0}=5, I recieved:

For E=0.05 and the same initial conditions, I recieved:

For E=0.075 and the same initial conditions, I recieved:

For E=0.1 and the same initial conditions, I recieved:

For E=0.15 and the same initial conditions, I recieved:

In all graphs elliptic fixed points can be found in the centers of the slliptic surfaces. chaotic regions can hardly be found. Therefore, I tried changing the initial xonditions.

For E=0.15 and larger y_{0}, I recieved:

This is still not chaotic.

For E=0.15 and larger p_{y0}, I recieved:

Finally we recieve a chaotic system.

3.

I choose to describe a problem discussed in an undergraduate course I have taken in the EE department - Introduction to Biological Signals and Systems.

In physiology, an action potential is a short-lasting event in which the electrical membrane potential of a cell rapidly rises and falls, following a consistent trajectory. The shape of a typical action potential is depicted in the following figure. The membrane potential remains near a baseline level until at some point in time it abruptly spikes upward and then rapidly falls.

The Hodgkinâ€“Huxley model, or conductance-based model, is a mathematical model that describes how action potentials in neurons are initiated and propagated, through a system of four first-order non-linear differential equations:

where V is the membrane potential, I is the externally applied current. V

We wish to solve the above set of equations for V(t), given a set of initial conditions and a given externally applied current.

I would chose to solve this problem using the fourth order Runge-Kutta method. This method is chosen because it is more accurate than the Euler method and since the equations are not linear, the implicit methods can be difficult. The steps for solving would be the following:

- Substitute the initial V(1) to calculate α(V(1)),β(V(1)) and substitute to get m(1),n(1),h(1).
- For each time step i: we shall define n
_{i}=dn/dt(i) and the same for m_{i}and h_{i}and calculate- k
_{1}=dt*HH(i,V(i),n(i),m(i),h(i)) where HH is the differential equation for dV/dt. - k
_{2}=dt*HH(i+dt/2,V(i)+k_{1}/2,n(i)+n_{1}/2,m(i)+m_{1}/2,h(i)+h_{1}/2) - k
_{3}=dt*HH(i+dt/2),V(i)+k_{2}/2,n(i)+n_{2}/2,m(i)+m_{2},h(i)+h_{2}) - k
_{4}=dt*HH(i+dt),V(i)+k_{3},n(i)+n_{3},m(i)+m_{3},h(i)+h_{3})

- k
- Create the next step for each variable:
- V(i+1)=V(i)+1/6*(k
_{1}+2*k_{2}+2*k_{3}+k_{4}) - n(i+1)=n(i)+1/6*(n
_{1}+2*n_{2}+2*n_{3}+n_{4}) - m(i+1)=m(i)+1/6*(m
_{1}+2*m_{2}+2*m_{3}+m_{4}) - h(i+1)=h(i)+1/6*(h
_{1}+2*h_{2}+2*h_{3}+h_{4})

- V(i+1)=V(i)+1/6*(k

For a given set of initial conditions and not too long a time vector, this should not be too long of a calculation, and I would plot V(t) at the end. But, if we wish to examine the response to changing external currents for a long period of time, I might choose to plot the progress of V(t) at every iteration.

The accuracy is the order of dt

Another problem might occur if the initial conditions and external current are such that will create chaos. Reduction of the above equations to two dimensions under certain assumptions can be examined in state-space representations and there isn't always a steady state. In real life this cannot happen (since n and h bring V down) but if the given conditions are not realistic one might receive chaos. This should not happen.