Solving system of equations with parameters

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!!

[Sophie Sophie]

I should say that I've written the following code and gotten an error:

Error using solve>getEqns (line 418)

List of equations must not be empty.

Error in solve (line 226)

[eqns,vars,options] = getEqns(varargin{:});

syms 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; 

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)';
tot_SERCA = '10000 = S0 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11'; 

[S0,S1,S2,S3,S4,S5,S6,S7,S8,S9,S10] = solve(dS0dt,dS1dt,dS2dt,dS3dt,dS4dt,dS5dt,dS6dt,dS7dt,dS8dt,dS9dt,dS10dt,dS11dt,tot_SERCA,'S0,S1,S2,S3,S4,S5,S6,S7,S8,S9,S10')

[Walter Roberson]

Which MATLAB release are you using? solve() cannot process strings since R2018a.

Answer 1

[madhan ravi]

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

[Sophie Sophie]

Thanks so much for the help!!

[madhan ravi]

Anytime :)

Answer 2

[Walter Roberson]

When you have more equations than unknowns then there is seldom a solution.

It looks like you want Ca to be a model parameter. However, there is no solution of the equations for arbitrary values of Ca, only for one specific value. Ca = -1.84683360434356e-07 -- Yes, it is negative.

syms S0 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11
syms Ca
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; 
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);
tot_SERCA = 10000 == S0 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11; 
sol = solve([dS0dt,dS1dt,dS2dt,dS3dt,dS4dt,dS5dt,dS6dt,dS7dt,dS8dt,dS9dt,dS10dt,dS11dt,tot_SERCA]);
disp(structfun(@double, sol, 'uniform', 0))

[Sophie Sophie]

Thanks Walter!

Can you clarify what youre doinng with Ca here and why I should be using it?

[Walter Roberson]

Your dS2dt involves Ca* k_S2_S3 where Ca has not been defined by you.

You define 12 variables explicitly S0 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 and you have 13 equations. By declaring Ca to be symbolic that makes it 13 equations in 13 unknowns, and since your equations are linear in the unknowns and your system is not singular, there is a unique solution.

MATLAB permits you to ask to solve for fewer variables than equations, but the answer is almost always empty.

[Walter Roberson]

If you leave out the dS2dt and tot_SERCA equation it is possible to solve the remaining set in terms of Ca and S11, the two variables you did not solve for. But then the tot_SERCA reduces to

10000 == (155526453553967425006438293929*S11)*1/4860228527020093611217586850

which forces S11 to a particular value, and the dS2dt reduces to

0 = -375077189670196373016559832188/20250952195917056713406611875*S11 - 235575680000/2349*S11*Ca

and once S11 has been forced, this forces the value for Ca as well.

[Sophie Sophie]

Ah!!!! This is 100% a typo.

Ca should have been Ca_cyt_conc and is therefore a known value. It has been corrected below. I should also have been solving for S11..

Therefore there are 13 equations and 12 unknowns

BTW, I am running 2016b

syms 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; 
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_cyt_conc * 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)';
tot_SERCA = '10000 = S0 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11'; 
[S0,S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11] = solve(dS0dt,dS1dt,dS2dt,dS3dt,dS4dt,dS5dt,dS6dt,dS7dt,dS8dt,dS9dt,dS10dt,dS11dt,tot_SERCA,'S0,S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11')

When I run this, I get the following..

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Support of character vectors that are not valid variable names or define

a number will be removed in a future release. To create symbolic expressions,

first create symbolic variables and then use operations on them.

> In sym>convertExpression (line 1559)

In sym>convertChar (line 1464)

In sym>tomupad (line 1216)

In sym (line 179)

In solve>getEqns (line 405)

In solve (line 225)

In SERCA_ODE_SS (line 43)

Warning: Do not specify equations and variables as character vectors. Instead,

create symbolic variables with syms.

> In solve>getEqns (line 445)

In solve (line 225)

In SERCA_ODE_SS (line 43)

S0 =

Empty sym: 0-by-1

S1 =

Empty sym: 0-by-1

S2 =

Empty sym: 0-by-1

S3 =

Empty sym: 0-by-1

S4 =

Empty sym: 0-by-1

S5 =

Empty sym: 0-by-1

S6 =

Empty sym: 0-by-1

S7 =

Empty sym: 0-by-1

S8 =

Empty sym: 0-by-1

S9 =

Empty sym: 0-by-1

S10 =

Empty sym: 0-by-1

S11 =

Empty sym: 0-by-1

[Walter Roberson]

With the revised equation, the solution for the first 12 equations is unique, and it is that all of the variables should be exactly 0. There is no solution in which the sum of them is anything non-zero.

If you have 13 linear equations in 12 variables, the only time there can be a solution is if one or more of the equations are redundant. (Nonlinear equations could potentially have solutions in such a situation though.)

[Sophie Sophie]

Does this mean that the solution to the state occupancies for S0, S1, .... S11 should all be zero?

I am super confused by how this could be. The model runs in a biochemical solver so I'm not certain how this is possible.

[Sophie Sophie]

Thanks so much for your patience and your help.

[Walter Roberson]

The solution of all 0 is correct provided that you have given the correct equations. You can work step by step.

eqn := [0 = 450*S1 - (1378*S0)*1/5 + 250*S11, 0 = (28*S0)*1/5 - 570*S1 + 25*S2,

0 = 120*S1 - 39*S2 + 16*S3, 0 = 14*S2 - 100016*S3 + 40*S4,

0 = 100000*S3 - 175*S4 + 765*S5, 0 = 135*S4 - 780*S5 + (27*S6)*1/10,

0 = 15*S5 - (21*S6)*1/5 + 2*S7,

0 = (3*S6)*1/2 - (6157266382182995*S7)*1/17592186044416 + 156*S8,

0 = 350*S7 - 316*S8 + 250*S9, 0 = 117*S10 - 390*S9,

0 = 300*S9 - 177*S10 + 60*S11, 0 = 270*S0 + 60*S10 - 310*S11];

[ 1378

eqn := [0 = 450 S1 - ---- S0 + 250 S11,

[ 5

28

0 = -- S0 - 570 S1 + 25 S2, 0 = 120 S1 - 39 S2 + 16 S3,

5

0 = 14 S2 - 100016 S3 + 40 S4, 0 = 100000 S3 - 175 S4 + 765 S5,

27 21

0 = 135 S4 - 780 S5 + -- S6, 0 = 15 S5 - -- S6 + 2 S7,

10 5

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

]

0 = 300 S9 - 177 S10 + 60 S11, 0 = 270 S0 + 60 S10 - 310 S11]

]

S0_ := solve(eqn[1], S0);

1125 625

S0_ := ---- S1 + --- S11

689 689

eqn2 := simplify(eval(eqn[2 .. -1], S0 = S0_));

[ 386430 3500

eqn2 := [0 = - ------ S1 + ---- S11 + 25 S2,

[ 689 689

0 = 120 S1 - 39 S2 + 16 S3, 0 = 14 S2 - 100016 S3 + 40 S4,

27

0 = 100000 S3 - 175 S4 + 765 S5, 0 = 135 S4 - 780 S5 + -- S6,

10

21

0 = 15 S5 - -- S6 + 2 S7,

5

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

303750 44840 ]

0 = ------ S1 - ----- S11 + 60 S10]

689 689 ]

S1_ := solve(eqn2[1], S1);

350 3445

S1_ := ----- S11 + ----- S2

38643 77286

eqn3 := simplify(eval(eqn2[2 .. -1], S1 = S1_));

[ 14000 433459

eqn3 := [0 = ----- S11 - ------ S2 + 16 S3,

[ 12881 12881

0 = 14 S2 - 100016 S3 + 40 S4, 0 = 100000 S3 - 175 S4 + 765 S5,

27 21

0 = 135 S4 - 780 S5 + -- S6, 0 = 15 S5 - -- S6 + 2 S7,

10 5

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

786860 253125 ]

0 = - ------ S11 + ------ S2 + 60 S10]

12881 12881 ]

S2_ := solve(eqn3[1], S2);

14000 206096

S2_ := ------ S11 + ------ S3

433459 433459

eqn4 := simplify(eval(eqn3[2 .. -1], S2 = S2_));

[ 196000 43349950000

eqn4 := [0 = ------ S11 - ----------- S3 + 40 S4,

[ 433459 433459

27

0 = 100000 S3 - 175 S4 + 765 S5, 0 = 135 S4 - 780 S5 + -- S6,

10

21

0 = 15 S5 - -- S6 + 2 S7,

5

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

26203540 4050000 ]

0 = - -------- S11 + ------- S3 + 60 S10]

433459 433459 ]

S3_ := solve(eqn4[1], S3);

14 433459

S3_ := ------- S11 + ---------- S4

3096425 1083748750

eqn5 := simplify(eval(eqn4[2 .. -1], S3 = S3_));

[ 56000 117048105

eqn5 := [0 = ------ S11 - --------- S4 + 765 S5,

[ 123857 866999

27 21

0 = 135 S4 - 780 S5 + -- S6, 0 = 15 S5 - -- S6 + 2 S7,

10 5

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

7487420 3240 ]

0 = - ------- S11 + ------ S4 + 60 S10]

123857 866999 ]

S4_ := solve(eqn5[1], S4);

78400 14738983

S4_ := -------- S11 + -------- S5

23409621 2601069

eqn6 := simplify(eval(eqn5[2 .. -1], S4 = S4_));

[ 392000 13023705 27

eqn6 := [0 = ------ S11 - -------- S5 + -- S6,

[ 867023 867023 10

21

0 = 15 S5 - -- S6 + 2 S7,

5

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

52413380 18360 ]

0 = - -------- S11 + ------ S5 + 60 S10]

867023 867023 ]

S5_ := solve(eqn6[1], S5);

78400 7803207

S5_ := ------- S11 + -------- S6

2604741 43412350

eqn7 := simplify(eval(eqn6[2 .. -1], S5 = S5_));

[ 392000 13056753

eqn7 := [0 = ------ S11 - -------- S6 + 2 S7,

[ 868247 8682470

3 6157266382182995

0 = - S6 - ---------------- S7 + 156 S8,

2 17592186044416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

52486820 16524 ]

0 = - -------- S11 + ------- S6 + 60 S10]

868247 4341235 ]

S6_ := solve(eqn7[1], S6);

3920000 17364940

S6_ := -------- S11 + -------- S7

13056753 13056753

eqn8 := simplify(eval(eqn7[2 .. -1], S6 = S6_));

[ 1960000 26645225141557261584225

eqn8 := [0 = ------- S11 - ----------------------- S7 + 156 S8,

[ 4352251 76565609303995580416

0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,

0 = 300 S9 - 177 S10 + 60 S11,

263095060 22032 ]

0 = - --------- S11 + ------- S7 + 60 S10]

4352251 4352251 ]

S7_ := solve(eqn8[1], S7);

1379227385882214400 3981411683807770181632

S7_ := ---------------------- S11 + ---------------------- S8

1065809005662290463369 8881741713852420528075

eqn9 := simplify(eval(eqn8[2 .. -1], S7 = S7_));

[ 482729585058775040000

eqn9 := [0 = ---------------------- S11

[ 1065809005662290463369

56525451689785812932020

- ----------------------- S8 + 250 S9, 0 = 117 S10 - 390 S9,

355269668554096821123

7158723507703201131260

0 = 300 S9 - 177 S10 + 60 S11, 0 = - ---------------------- S11

118423222851365607041

6718244744129937408 ]

+ ---------------------- S8 + 60 S10]

2960580571284140176025 ]

S8 := 'S8';

S8_ := solve(eqn9[1], S8);

S8 := S8

24136479252938752000 8881741713852420528075

S8_ := ---------------------- S11 + ---------------------- S9

8478817753467871939803 5652545168978581293202

eqn10 := simplify(eval(eqn9[2 .. -1], S8 = S8_));

[

eqn10 := [0 = 117 S10 - 390 S9, 0 = 300 S9 - 177 S10 + 60 S11, 0 =

[

170849106624052570889180 10077367116194906112

- ------------------------ S11 + ---------------------- S9

2826272584489290646601 2826272584489290646601

]

+ 60 S10]

]

S9_ := solve(eqn10[1], S9);

3

S9_ := -- S10

10

eqn11 := simplify(eval(eqn10[2 .. -1], S9 = S9_));

[

eqn11 := [0 = -87 S10 + 60 S11, 0 =

[

170849106624052570889180 847896891397461486339468 ]

- ------------------------ S11 + ------------------------ S10]

2826272584489290646601 14131362922446453233005 ]

S10_ := solve(eqn11[1], S10);

20

S10_ := -- S11

29

eqn12 := simplify(eval(eqn11[2 .. -1], S10 = S10_));

[ 1563036526507678610428348 ]

eqn12 := [0 = - ------------------------- S11]

[ 81961904950189428751429 ]

S11_ := solve(eqn12[1], S11);

S11_ := 0

And now with S11 being 0, S10 is forced to be 0, and everything else back substitutes to 0. No escaping it.

I would have a question for you, though: the names of your variables suggest that you have differential equations rather than simple linear equations. Are you sure you should be using plain solve() and not dsolve() of something or other?