% ===============================================================
% == Multi Affine Minimal Model Simulation Module              ==
% ==                                                           ==
% ==   Supplement of the work:                                 ==
% ==                                                           ==
% ==     From Cardiac Cells to Genetic Regulatory Networks     ==
% ==     Parameter Estimation for Action Potential Duration    ==
% ==                                                           ==
% ==   Author:                                                ==
% ==                                                           ==
% ==     E. Bartocci                                           ==
% ==                                                           == 
% ==   Date:  11/05/10                                         ==
% ==                                                           ==
% ==   Free distribution with authors permission               ==
% ==                                                           ==
% ==   SUNY Stony Brook, Stony Brook, NY                       ==
% ==                                                           == 
% ===============================================================  
function [out] = sim_cell_model ()

         global u_partition;
         global v_partition;
         global w_partition;
         global s_partition;
         global q_partition;
         global dt;
         
         u_partition = [0 0.00322522522522523 0.0059 0.006 0.0197777777777778 0.0293353353353353 0.0520500500500501 0.0618558558558559 0.0765025025025025 0.129 0.13 0.176116116116116 0.219509509509509 0.260690690690691 0.29 0.3 0.426376376376376 0.541491491491492 0.644094094094094 0.734184184184184 0.810510510510511 0.874324324324324 0.925625625625626 0.976926926926927 1.03948948948949  1.11581581581582  1.20590590590591 1.30850850850851 1.42362362362362 1.55];
         v_partition = [0 0.95 1];
         w_partition = [0 0.95 1];
         s_partition = [0 0.01 1];
         q_partition = [1.0  1.021426  1.02251  600.0];

         ut1     = [];
         vt1     = [];
         wt1     = [];
         st1     = [];
         stimt1  = [];
      
         u = 0.0;
         v = 1.0;
         w = 1.0;
         s = 0.0;
         q = 0.0;
         dt = 0.02;

         for i=1:30000   
             [u,v,w,s,q, stim] = compute_PMA_step (u,v,w,s,q);
             ut1(i) =    u;
             vt1(i) =    v;
             wt1(i) =    w;
             st1(i) =    s;
             stimt1(i)=  stim;     
             
         end
         
         %hold off;
         plot(linspace(0,600, 30000), [ut1; vt1; wt1; st1]);
         xlabel('Time in milliseconds (ms)');
         %plot(linspace(0,30, 30000), ut3);
         
 end

function [u,v,w,s,q, stim] = compute_PMA_step (u,v,w,s,q)
         global u_partition;
         global v_partition;
         global w_partition;
         global s_partition;
         global q_partition;      
         global dt;

         prod1 = 0.0166666666666667   * R (u, [1 3 4 30],       u_partition, [1.0 1.0 0.0 0.0]) * (1 - v); 
         degr1 = 0.000869565217391304 * R (u, [1 3 4 15 16 30], u_partition, [0.0 0.0 1.0 1.0 0.0 0.0]) * v;
         degr2 = 0.689369915896870    * R (u, [1 15 16 30],     u_partition, [0.0 0.0 1.0 1.0]) * v;
         
         v = v + (prod1 - degr1 - degr2)*dt;
      
         prod1 = 0.005 * R(u, [1 3 4 30],       u_partition, [1.0 1.0 0.0 0.0]) * (1 - w);
         prod2 = 0.005 * R(u, [1 3 4 10 11 30], u_partition, [0.0 0.0 0.94 0.94 0.0 0.0]);
         prod3 = -0.071428571428571 * R(u, [1 3 4 30], u_partition, [0.0 0.005999999 0.0 0.0]);
         degr1 = 0.005 * R(u, [1 3 4 10 11 30], u_partition, [0.0 0.0 1.0 1.0 0.0 0.0]) * w;
         degr2 = 0.0025 * R(u, [1 10 11 30], u_partition, [0.0 0.0 1.0 1.0]) * w;
            
         w = w + (prod1 + prod2 + prod3 - degr1 - degr2) * dt;
      
         prod1 = 0.365737692926633  * R (u, [1 2 3 4 5 6 7 8 9 10 11 30], u_partition, [ 0.0215530430802808   0.0218404832501762	 0.0220907746541382   0.0220907747448439   0.0233756631454813	0.0243095439036467	 0.0266773222136253	  0.0277674867419459   0.0294768269421839   0.0366287295501185   0.0	0.0 ]);
         prod2 = 0.0625000000000000 * R (u, [1 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30], u_partition, [ 0.0   0.0	        0.0366287443664546   0.0441092881834653   0.0524622441859381   0.0617535327747815	0.0720368814519242   0.0720371621320020   0.116584384855182   0.176268665677548   0.247679421847480   0.324590012929517   0.398366018652226   0.463978368245852	  0.517759355285649	  0.571132010712343	  0.633940429340140	  0.704670344322184	  0.776937543225223	  0.842734428913138	  0.896790494026689	  0.936593943204126]);
         degr1 = 0.365737692926633  * R (u, [1 10 11 30], u_partition, [1.0 1.0 0.0 0.0 ] ) * s;
         degr2 = 0.0625000000000000 * R (u, [1 10 11 30], u_partition, [0.0 0.0 1.0 1.0 ]) * s;
            
         s = s + (prod1 + prod2 - degr1 - degr2) * dt;
      
         q = q + dt;
      
         
         %stim = 1; %R (q, [1  2  3  4 ], q_partition, [1  1  0  0 ]); 
      
         %This line is just for having a small stimulus instead of the
         %original model
         stim = R (q, [1  2  3  4 ], q_partition, [1  1  0  0 ]); 
         jsi= 0.529801324503311 * R (u, [1 10 11 30], u_partition, [0  0  1  1])* R (w, [1 3], w_partition, [0  1]) * R (s, [1 3], s_partition, [0  1]);
         
         jso1 = 0.0025 * R (u, [1 2 3 4 30], u_partition, [0.0 0.00322522522522523 0.005999999 0.0 0.0]);
         jso2 = 1.66666667 * R (u, [1 3 4 5 6 7 8 9 10 11 30], u_partition, [0.0 0.0 0.006 0.0197777777777778 0.0293353353353353 0.0520500500500501 0.0618558558558559 0.0765025025025025 0.1299999  0.0  0.0]);
         jso3 = R(u, [ 1 10 11 12 13 14 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30], u_partition, [0.0 0.0 0.0371347486842788  0.0379245744529007	  0.0388154159747791  0.0398184365529367	0.0409446382840491 0.0460443086347673  0.0535437875115028  0.0637664813337137  0.0767333519436878  0.0916925245941391  0.107774773920785  0.123549833090718  0.142260797964051	0.169583455745153  0.210423641565343  0.270114601421950  0.352978734119542	0.460653839673183  0.585471819608690]);
         
         jfi= R(u, [1 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30], u_partition, [0.0 0.0	1.29090438149497	2.21405658092345	2.83378974570166	3.22013022213222	3.43197415269487	3.52779068995285	3.55113280538706	3.52662358418114	3.43197415269487	3.22013022213222	2.83378974570166	2.21405658092345	1.29090438149497	0.0]) * R (v, [1 3], v_partition, [0  1]);
       
         u = u  + (jfi-jso1-jso2-jso3+jsi+stim)*dt;
end


function [y] = R (x, ixs, xs, ys)
     for xi=2:size(ixs,2)
             if x < xs(ixs(xi))
                break; 
             end
     end
         
     x1 = xs(ixs(xi-1));
     x2 = xs(ixs(xi));
     y1 = ys(xi-1);
     y2 = ys(xi);
     y = y1 + ((y2 - y1) / (x2 - x1)) * (x -x1);
end

function [y] = H (x, theta, y1, y2)
     if (x - theta) < 0
         y = y2;
     else
         y = y1; 
     end
end

function [y] = S (x, k, theta, a, b)
     y = a + (b - a) * 0.5 * (1.+tanh(k * (x- theta)));
end

