% =============================================================== % == 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