MATLAB: Solving system of equations with parameters

parameterssystem of equations

Hi MATLAB community,
I am hoping someone can point me towards some code that I can use to solve a system of 13 equations with 12 unknowns.
All my differential equations are equal to 0 except for the last equation which is the sum of all my states.
State 0 is named S0 etc. Theres a total of 12 states.
I am solving for the values S0, S1, S2, S3, S4 , S5, S6, S7, S8. S9. S10, and S11
Here are the set of 13 equations with 12 unknowns and the model parameters below:
dS0dt  = 0 = ((k_S11_S0 * S11) + (k_S1_S0 * S1)) - ((Ca_cyt_conc* k_S0_S1 + Pi_conc * k_S0_S11) * S0);
dS1dt  = 0 = ((Ca_cyt_conc* k_S0_S1 * S0) + (k_S2_S1 * S2)) - ((k_S1_S0 + k_S1_S2) * S1);
dS2dt  = 0 = ((k_S1_S2 * S1) + (k_S3_S2 * S3)) - ((k_S2_S1 + Ca* k_S2_S3) * S2);
dS3dt  = 0 = ((Ca_cyt_conc* k_S2_S3 * S2) + (k_S4_S3 * S4)) - ((k_S3_S2 + MgATP_conc * k_S3_S4) * S3);
dS4dt  = 0 = ((MgATP_conc * k_S3_S4 * S3) + (k_S5_S4 * S5)) - ((k_S4_S3 + k_S4_S5) * S4);
dS5dt  = 0 = ((k_S4_S5 * S4) + (k_S6_S5 * MgADP_conc * S6)) - ((k_S5_S4 + k_S5_S6) * S5);
dS6dt  = 0 = ((k_S5_S6 * S5) + (k_S7_S6 * S7)) - ((k_S6_S5 * MgADP_conc + k_S6_S7) * S6);
dS7dt  = 0 = ((k_S6_S7 * S6) + (Ca_sr_conc * k_S8_S7 * S8)) - ((k_S7_S6 * MgADP_conc + k_S7_S8) * S7);
dS8dt  = 0 = ((k_S7_S8 * S7) + (k_S9_S8 * S9)) - ((Ca_sr_conc * k_S8_S7 + k_S8_S9) * S8);
dS9dt  = 0 = ((k_S8_S9 * S9) + (Ca_sr_conc * k_S10_S9 * S10)) - ((k_S9_S8 + k_S9_S10) * S9);
dS10dt = 0 = ((k_S9_S10 * S9) + (k_S11_S10 * S11)) - ((k_S10_S11 + Ca_sr_conc * k_S10_S9) * S10);
dS11dt = 0 = ((k_S10_S11 * S10) + (Pi_conc * k_S0_S11 * S0)) - ((k_S11_S0 + k_S11_S10) * S11);
10000 = S0 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11; 
Ca_sr_conc  = 1.3e-3;
MgATP_conc  = 5e-3;  
MgADP_conc  = 36e-6; 
Pi_conc     = 1e-3;  
Ca_cyt_conc = 140e-9; 
k_S0_S1   = 4e7;   
k_S1_S0   = 450;   
k_S1_S2   = 120;   
k_S2_S1   = 25;    
k_S2_S3   = 1e8;   
k_S3_S2   = 16;    
k_S3_S4   = 2e7;  
k_S4_S3   = 40;   
k_S4_S5   = 135;   
k_S5_S4   = 765;   
k_S5_S6   = 15;    
k_S6_S5   = 7.5e4;
k_S6_S7   = 1.5; 
k_S7_S6   = 2.0;
k_S7_S8   = 350; 
k_S8_S7   = 1.2e5;
k_S8_S9   = 160;   
k_S9_S8   = 250;   
k_S9_S10  = 300;
k_S10_S9  = 9e4;   
k_S10_S11 = 60;    
k_S11_S10 = 60;    
k_S11_S0  = 250;   
k_S0_S11  = 2.7e5;
Thanks in advance for your help!!

Best Answer

Like sir Walter says I think you should be using ode solvers to solve the problem because it clearly is a system of odes:
tspan=[0 1];  % time span 
ic=[0;1;0;0;0;0;0;0;0;0;0;0]; % initial conditions
[t,x]=ode45(@myod,tspan,ic);  % function call
plot(t,x(:,1)) % likewise another 11 solutions can be plotted
function dSdt = myod(t,x)  % function definition
 S0=x(1);
 S1=x(2);
 S2=x(3);
 S3=x(4);
 S4=x(5);
 S5=x(6);
 S6=x(7);
 S7=x(8);
 S8=x(9);
 S9=x(10);
 S10=x(11);
 S11=x(12);
Ca_sr_conc  = 1.3e-3;
MgATP_conc  = 5e-3;  
MgADP_conc  = 36e-6; 
Pi_conc     = 1e-3;  
Ca=1e-3;  % change it to your value
Ca_cyt_conc = 140e-9; 
k_S0_S1   = 4e7;   
k_S1_S0   = 450;   
k_S1_S2   = 120;   
k_S2_S1   = 25;    
k_S2_S3   = 1e8;   
k_S3_S2   = 16;    
k_S3_S4   = 2e7;  
k_S4_S3   = 40;   
k_S4_S5   = 135;   
k_S5_S4   = 765;   
k_S5_S6   = 15;    
k_S6_S5   = 7.5e4;
k_S6_S7   = 1.5; k_S7_S6   = 2.0;
k_S7_S8   = 350; 
k_S8_S7   = 1.2e5;
k_S8_S9   = 160;   
k_S9_S8   = 250;   
k_S9_S10  = 300;
k_S10_S9  = 9e4;   
k_S10_S11 = 60;    
k_S11_S10 = 60;    
k_S11_S0  = 250;   
k_S0_S11  = 2.7e5; 
dSdt  = [ ((k_S11_S0 * S11) + (k_S1_S0 * S1)) - ((Ca_cyt_conc* k_S0_S1 + Pi_conc * k_S0_S11) * S0);
((Ca_cyt_conc* k_S0_S1 * S0) + (k_S2_S1 * S2)) - ((k_S1_S0 + k_S1_S2) * S1);
((k_S1_S2 * S1) + (k_S3_S2 * S3)) - ((k_S2_S1 + Ca* k_S2_S3) * S2);
((Ca_cyt_conc* k_S2_S3 * S2) + (k_S4_S3 * S4)) - ((k_S3_S2 + MgATP_conc * k_S3_S4) * S3);
((MgATP_conc * k_S3_S4 * S3) + (k_S5_S4 * S5)) - ((k_S4_S3 + k_S4_S5) * S4);
((k_S4_S5 * S4) + (k_S6_S5 * MgADP_conc * S6)) - ((k_S5_S4 + k_S5_S6) * S5);
((k_S5_S6 * S5) + (k_S7_S6 * S7)) - ((k_S6_S5 * MgADP_conc + k_S6_S7) * S6);
((k_S6_S7 * S6) + (Ca_sr_conc * k_S8_S7 * S8)) - ((k_S7_S6 * MgADP_conc + k_S7_S8) * S7);
((k_S7_S8 * S7) + (k_S9_S8 * S9)) - ((Ca_sr_conc * k_S8_S7 + k_S8_S9) * S8);
((k_S8_S9 * S9) + (Ca_sr_conc * k_S10_S9 * S10)) - ((k_S9_S8 + k_S9_S10) * S9);
((k_S9_S10 * S9) + (k_S11_S10 * S11)) - ((k_S10_S11 + Ca_sr_conc * k_S10_S9) * S10);
((k_S10_S11 * S10) + (Pi_conc * k_S0_S11 * S0)) - ((k_S11_S0 + k_S11_S10) * S11)];
end
Related Question