Chapter 3 Using Neural Networks to Model Biological Objects:

3.1 A Neural Network to Model An Action Pulse:

We set up a simple 1 - 7 - 7 - 1 network using sigmoidal transfer functions on the input and hidden layers and a linear transfer function on the single output neuron. The plan is to do this: first, generate the action pulse, then find the portion of the action pulse we want to model and then do the modeling. Here we use a new updated HH solver given below:

SolveSimpleHH2.m

 1 function [t,I_NA,I_K,g_NA,g_K,V,m_NA,h_NA,n_K] = ...
 2           SolveSimpleHH2(fname,t0,t1,y0,T,...
 3     Kconc,g_Kmax,...
 4     NAconc,g_NAmax)
 5 
 6 %
 7 % We use a simple Runge-Kutta Fehlberg scheme
 8 % using the Matlab function below:
 9 %
10 %function [t,y] = ode45('fname',tspan,y0,options)
11 %
12 %        Gives approximate solution to 
13 %          y'(t) = f(t,y(t))
14 %          y(tmin) = y0
15 %        using a 4-5th order RKF method 
16 %
17 % t0     initial time
18 % t1     final time
19 % y0     initial state
20 % tspan  [t0 t1]--column vector
21 %
22 % t      time values returned by ODE45 algorithm
23 %        --not necessarily uniformly spaced
24 %
25 % T      = Temperature in Celsius
26 % NAconc = [Na_in,NA_out]--column vector 
27 % g_Na   = sodium conductance
28 % m_NA   = activation for sodium
29 % h_NA   = inactivation for sodium
30 %
31 % Kconc  = [Kin,Kout]--column vector 
32 % g_K  = potassium conductance
33 % n_K  = activation for potassium  
34 % V    = membrane voltage 
35 %
36 tspan = [t0;t1];
37 options = odeset('RelTol',1.0e-6);
38 [t,y] = ode45(fname,tspan,y0,options);
39 %
40 % Set solution variables
41 %
42 V    = y(:,1); 
43 m_NA = y(:,2);
44 h_NA = y(:,3);
45 n_K  = y(:,4);
46 %
47 %  Compute g_Na and G_K
48 %
49 g_NA = g_NAmax*m_NA.*m_NA.*m_NA.*h_NA;
50 g_K  = g_Kmax*n_K.*n_K.*n_K.*n_K;
51 E_NA = Nernst(1,T,NAconc(1),NAconc(2));
52 I_NA = g_NA.*(V - E_NA);
53 E_K = Nernst(1,T,Kconc(1),Kconc(2));
54 I_K = g_NA.*(V - E_K);

syntax highlighted by Code2HTML, v. 0.8.8b

and the actual session is as follows:

runtime3.prompts

  1 
  2                               < M A T L A B >
  3                   Copyright 1984-2001 The MathWorks, Inc.
  4                        Version 6.1.0.450 Release 12.1
  5                                 May 18 2001
  6 
  7  
  8   To get started, select "MATLAB Help" from the Help menu.
  9  
 10 >> path(path,'/local/petersj/BioInfo/HH');
 11 >> path(path,'/local/petersj/BioInfo/Nernst');
 12 >> [t,I_NA,I_K,g_NA,g_K,V,m_NA,h_NA,n_K] = ...
 13 SolveSimpleHH2('simpleHH2',0,128,[-70.0;0.2;0.8;0.2],9.3,[400.0;20.11],...
 14 36.0,[50.0;491.0],120.0);
 15 >> plot(t,V,'g-');
 16 %
 17 %  Just do 20 seconds
 18 %
 19 >> [t,I_NA,I_K,g_NA,g_K,V,m_NA,h_NA,n_K] = ...
 20 SolveSimpleHH2('simpleHH2',0,20,[-70.0;0.2;0.8;0.2],9.3,[400.0;20.11],...
 21 36.0,[50.0;491.0],120.0);
 22 >> plot(t,V,'g-');
 23 %
 24 % Process data
 25 %
 26 >> length(t)
 27 
 28 ans =
 29 
 30    889
 31 
 32 >> length(V)
 33 
 34 ans =
 35 
 36    889
 37 
 38 >> plot(t(45:600),V(45:600),'g-');
 39 >> plot(t(445:600),V(445:600),'g-');
 40 >> plot(t(345:600),V(345:600),'g-');
 41 >> plot(t(245:600),V(245:600),'g-');
 42 >> plot(t(205:600),V(205:600),'g-');
 43 >> plot(t(185:600),V(185:600),'g-');
 44 >> plot(t(165:650),V(165:650),'g-');
 45 >> plot(t(155:690),V(155:690),'g-');
 46 >> plot(t(155:790),V(155:790),'g-');
 47 >> plot(t(135:820),V(135:820),'g-');
 48 >> plot(t(135:840),V(135:840),'g-');
 49 >> plot(t(135:889),V(135:889),'g-');
 50 %
 51 % OK found pulse
 52 %
 53 %
 54 % Now find NN model
 55 %
 56 >> net = newff([tmin tmax], [1 7 7 1],{'tansig','tansig','tansig','purelin'});
 57 >> Y = sim(net,P);
 58 >> plot(P,T,'*',P,Y);
 59 >> net.trainParam.epochs = 400;
 60 >> net = train(net,P,T);
 61 TRAINLM, Epoch 0/400, MSE 2172.88/0, Gradient 439928/1e-10
 62 TRAINLM, Epoch 7/400, MSE 2053.53/0, Gradient 1.47937e-09/1e-10
 63 TRAINLM, Maximum MU reached, performance goal was not met.
 64 %
 65 % This trains badly -- note that the outpu range is huge
 66 % so we have a real hard time hitting targets.
 67 %
 68 >> Min = min(T)
 69 
 70 Min =
 71 
 72   -72.5083
 73 
 74 >> Max = max(T)
 75 
 76 Max =
 77 
 78    50.6548
 79 
 80 >> D = Max-Min
 81 
 82 D =
 83 
 84   123.1631
 85 
 86 >> S = Min+Max
 87 
 88 S =
 89 
 90   -21.8535
 91 
 92 >> T2 = (2/D)*(T-0.5*S);
 93 >> min(T2)
 94 
 95 ans =
 96 
 97     -1
 98 
 99 >> max(T2)
100 
101 ans =
102 
103      1
104 
105 >>
106  >> net2 = newff([tmin tmax],[1 7 7 1],{'tansig','tansig','tansig','purelin'});
107 >> Y = sim(net2,P);
108 >> plot(P,T2,'*',P,Y,'g-');
109 >> plot(P,T2,'*',P,Y,'g-');
110 >> net2.trainParam.epochs = 400;
111 >> net2 = train(net2,P,T2);
112 TRAINLM, Epoch 0/400, MSE 0.323757/0, Gradient 5621.48/1e-10
113 TRAINLM, Epoch 25/400, MSE 0.000560481/0, Gradient 3.8885/1e-10
114 TRAINLM, Epoch 50/400, MSE 0.000397259/0, Gradient 1.0903/1e-10
115 TRAINLM, Epoch 75/400, MSE 0.000332295/0, Gradient 2.9831/1e-10
116 TRAINLM, Epoch 100/400, MSE 0.000321657/0, Gradient 4.91938/1e-10
117 TRAINLM, Epoch 125/400, MSE 0.000312525/0, Gradient 0.878092/1e-10
118 TRAINLM, Epoch 150/400, MSE 0.000165219/0, Gradient 4.38306/1e-10
119 TRAINLM, Epoch 175/400, MSE 1.57439e-05/0, Gradient 4.13758/1e-10
120 TRAINLM, Epoch 200/400, MSE 6.21834e-06/0, Gradient 0.416302/1e-10
121 TRAINLM, Epoch 225/400, MSE 5.67268e-06/0, Gradient 1.48881/1e-10
122 TRAINLM, Epoch 250/400, MSE 4.50269e-06/0, Gradient 0.851653/1e-10
123 TRAINLM, Epoch 275/400, MSE 3.32058e-06/0, Gradient 0.573192/1e-10
124 TRAINLM, Epoch 300/400, MSE 2.6245e-06/0, Gradient 0.0213724/1e-10
125 TRAINLM, Epoch 325/400, MSE 1.92723e-06/0, Gradient 0.827397/1e-10
126 TRAINLM, Epoch 350/400, MSE 1.6671e-06/0, Gradient 0.0257573/1e-10
127 TRAINLM, Epoch 375/400, MSE 1.61861e-06/0, Gradient 0.042919/1e-10
128 TRAINLM, Epoch 400/400, MSE 1.5497e-06/0, Gradient 0.492028/1e-10
129 TRAINLM, Maximum epoch reached, performance goal was not met.
130 
131 >> Y = sim(net2,P);
132 >> plot(P,T2,'*',P,Y,'g-');
133 >> Q = linspace(0,40,1000);
134 >> Z = sim(net2,Q);
135 >> plot(P,T2,'*',Q,Z,'g-');

syntax highlighted by Code2HTML, v. 0.8.8b

3.1.1 Prior To Training:

Before updating the parameters via back propagation to meet the desired targets, we see the network performance is poor as seen in Figure 3.1.1:

Figure 3.1: Prior To Training:

3.1.2 Neural Network Output versus The Targets:

After training, we see the results in Figure 3.1.2:

Figure 3.2: After Training:

3.1.3 Neural Network Output versus The Targets: Testing

Here in Figuree 3.1.3, we see how the network performs on data it has not been trained on:

Figure 3.3: Testing:

3.2 Neural Models of Gates:

We will solve the standard HH model using neural networks for some of our computations.

We will model each of the a and b functions used in the HH gate model with a neural network.
We can then replace the a and b functions in the HH ODE function with the neural networks from above. This will require that you write new MatLab functions to handle the HH code.
Using your new neural net code, we can then compute the action potential as before. We can then do parametric studies on the action potential output using various parameters.

3.2.1 Using Neural Models in the Axon Pulse Generation:

Now we create a neural model for one of the a functions that we use in the computation of the axonal pulse. We will use a very primitive neural net creation function just to get started. We will model a_mNA from the simpleHH2 function; we do this with the CreateNN function below:

CreateNN.m

 1 function net = CreateNN(name,x)
 2 %
 3 % function = name of function to model
 4 % x = input data
 5 %
 6 n = length(x);
 7 y = zeros(1,n);
 8 for i=1:length(x)
 9   y(i) = feval(name,x(i));
10 end
11 low = min(x)
12 high = max(x)
13 
14 %
15 % Scale targets
16 %
17 %L = min(y)
18 %H = max(y)
19 %D = (H-L)/2.0;
20 %A = (H+L)/2.0;
21 %y2 = (y-A)/D;
22 %min(y2)
23 %max(y2) 
24 net = newff([low high],[ 1 7 7 1], {'tansig','tansig','tansig','purelin'});
25 Y = sim(net,x);
26 plot(x,y,'*',x,Y,'g-');
27 net.trainParam.epochs = 200;
28 net = train(net,x,y);
29 Y = sim(net,x);
30 plot(x,y,'*',x,Y,'g-');

syntax highlighted by Code2HTML, v. 0.8.8b This will create a simple neural model of our a. We will make this a global variable with the line global nnalpha_mNA; in our session as you will see. We then declare this variable to be global in our new simpleHHNN function in the same way. In the new simpleHHNN.m file below, note we have really only removed the old alpha_mNA declarations and replaced with the line alpha_mNA = sim(nnalpha_mNA,y(1)); Unfortunately, this makes the code run a lot slower. When we run our old simulation code we generate the following pulse: as seen in Figure ( 3.2.1):

Figure 3.4:

When we use the neural model, we get the same performanceas seen in Figure ( 3.2.1);

Figure 3.5:

Here is the log of our MatLab session so you can follow along yourself:

session120501.prompts

  1 %
  2 % Now let's model alpha_mNA
  3 %
  4 >> V = linspace(-100.0,100.0,200);
  5 
  6 %
  7 % Create global variable;
  8 %
  9 >> global nnalpha_mNA;
 10 
 11 %
 12 % We wrote a simple nn creation function to
 13 % see if we can do this
 14 % this will create a 1-7-7-1 NN with
 15 % 'tansig' neurons on all but the output which is
 16 % 'purelin'
 17 %
 18 
 19 >> nnalpha_mNA = CreateNN('alpha_mNA',V);
 20 
 21 %
 22 % echo min and max of input data
 23 %
 24 low =
 25 
 26   -100
 27 
 28 
 29 high =
 30 
 31    100
 32 
 33 
 34 TRAINLM, Epoch 0/200, MSE 27.4252/0, Gradient 12720.1/1e-10
 35 TRAINLM, Epoch 25/200, MSE 0.0395543/0, Gradient 183.413/1e-10
 36 TRAINLM, Epoch 50/200, MSE 0.0212297/0, Gradient 65.0449/1e-10
 37 TRAINLM, Epoch 75/200, MSE 0.0133573/0, Gradient 16.9445/1e-10
 38 TRAINLM, Epoch 100/200, MSE 0.00561581/0, Gradient 730.551/1e-10
 39 TRAINLM, Epoch 125/200, MSE 6.65846e-05/0, Gradient 35.1209/1e-10
 40 TRAINLM, Epoch 150/200, MSE 1.55648e-05/0, Gradient 16.9909/1e-10
 41 TRAINLM, Epoch 175/200, MSE 2.83548e-06/0, Gradient 0.607174/1e-10
 42 TRAINLM, Epoch 200/200, MSE 1.04892e-06/0, Gradient 3.14067/1e-10
 43 TRAINLM, Maximum epoch reached, performance goal was not met.
 44 
 45 %
 46 % so now we have created the neural model of alpha_mNA called
 47 % nnalpha_mNA which is global
 48 %
 49 
 50 %
 51 % Let's see if our solver routines still work
 52 %
 53 
 54 >> [t,I_NA,I_K,g_NA,g_K,V,m_NA,h_NA,n_K] = ...
 55 SolveSimpleHH2('simpleHH2',0,128,[-70.0;0.2;0.8;0.2],9.3,[400.0;20.11],...
 56 36.0,[50.0;491.0],120.0);
 57 
 58 %
 59 % lot's of division by zero things
 60 %
 61 
 62 Warning: Divide by zero.
 63 > In /local/petersj/BioInfo/HH/simpleHH2.m at line 92
 64   In /usr/commercial/matlab12.1/toolbox/matlab/funfun/ode45.m at line 311
 65   In /local/petersj/BioInfo/HH/SolveSimpleHH2.m at line 38
 66 Warning: Divide by zero.
 67 > In /local/petersj/BioInfo/HH/simpleHH2.m at line 115
 68   In /usr/commercial/matlab12.1/toolbox/matlab/funfun/ode45.m at line 311
 69   In /local/petersj/BioInfo/HH/SolveSimpleHH2.m at line 38
 70 Warning: Divide by zero.
 71 > In /local/petersj/BioInfo/HH/simpleHH2.m at line 92
 72   In /usr/commercial/matlab12.1/toolbox/matlab/funfun/ode45.m at line 320
 73   In /local/petersj/BioInfo/HH/SolveSimpleHH2.m at line 38
 74 Warning: Divide by zero.
 75 > In /local/petersj/BioInfo/HH/simpleHH2.m at line 115
 76   In /usr/commercial/matlab12.1/toolbox/matlab/funfun/ode45.m at line 320
 77   In /local/petersj/BioInfo/HH/SolveSimpleHH2.m at line 38
 78 Warning: Divide by zero.
 79 > In /local/petersj/BioInfo/HH/simpleHH2.m at line 141
 80   In /usr/commercial/matlab12.1/toolbox/matlab/funfun/ode45.m at line 320
 81   In /local/petersj/BioInfo/HH/SolveSimpleHH2.m at line 38
 82   
 83 %
 84 % Fix this by changing 'simpleHH2' to use reltol = 1.0e-9.
 85 %
 86 
 87 >> [t,I_NA,I_K,g_NA,g_K,V,m_NA,h_NA,n_K] = ...
 88 SolveSimpleHH2('simpleHH2',0,128,[-70.0;0.2;0.8;0.2],9.3,[400.0;20.11],...
 89 36.0,[50.0;491.0],120.0);
 90 
 91 %
 92 % this is the usual plot of the pulse
 93 %
 94 >> plot(t,V,'g-');
 95 
 96 %
 97 % now we generate the plot using the neural model of alpha_mNA
 98 %
 99 
100 >> [t,I_NA,I_K,g_NA,g_K,V,m_NA,h_NA,n_K] = ...
101 SolveSimpleHH2('simpleHHNN',0,128,[-70.0;0.2;0.8;0.2],9.3,[400.0;20.11],...
102 36.0,[50.0;491.0],120.0);
103 
104 %
105 % This took a LOT longer to compute because of all those
106 % neural model calculations.
107 % The plots look the same
108 %
109 
110 >> plot(t,V,'g-')

syntax highlighted by Code2HTML, v. 0.8.8b