How to solve 5 quadratic equations simultaneously

I am wondering if anyone has experience with using Matlab to solve 5 quadratic equations simultaneously – we are trying to model competitive binding with 5 proteins and 1 ligand. The program seems to work fine with up to 4 equations but crashes with 5 equations.

My codes are as follows:

(for 4)

syms LIDtot I1tot kdI1 x I2tot kdI2 y L3tot kdL3 z L4tot kdL4 w
S = solve ([(I1tot-x)*(LIDtot-x-y-z-w)==kdI1*x,(I2tot-y)*(LIDtot-x-y-z-w)==kdI2*y,(L3tot-z)*(LIDtot-x-y-z)==kdL3*z,(L4tot-w)*(LIDtot-x-y-z-w)==kdL4*w], [x,y,z,w])

This works fine.

(for 5)

syms LIDtot I1tot kdI1 x I2tot kdI2 y L3tot kdL3 z L4tot kdL4 w M4tot kdM4 v
S = solve ([(I1tot-x)*(LIDtot-x-y-z-w-v)==kdI1*x,(I2tot-y)*(LIDtot-x-y-z-w-v)==kdI2*y,(L3tot-z)*(LIDtot-x-y-z-w-v)==kdL3*z,(L4tot-w)*(LIDtot-x-y-z-w-v)==kdL4*w,(M4tot-v)*(LIDtot-x-y-z-w-v)==kdM4*v], [x,y,z,w,v])

This crashes.

Any help and advice would be much appreciated. Thanks.

Answer 1

[Roger Stafford]

For that particular set of equations you could take steps to make things easier. Introduce a new unknown 't':

t = LIDtot-x-y-z-w-v;

Then each of your five equations can easily be solved for the corresponding x, y, z, w, or v in terms of t. Finally, introducing those expressions into

t = LIDtot-x-y-z-w-v;

will give you a single equation in 't'. I think this equation will simplify to a fifth degree polynomial equation in t, so there may not exist an explicit solution for it, but it presents a direct way of numerically solving your equations using matlab's 'roots' functions in combination with the above expressions for x, y, z, w, and u, for any given numerical values of LIDtot, I1tot, kdI1, etc.

Also this perhaps shows why you succeeded with only four equations and four unknowns since that would have resulted in a quartic equation in t which possesses an explicit solution.

[Roger Stafford]

As Walter has pointed out, strictly speaking, the single equation in t will simplify to a sixth degree polynomial, but t is a factor in all its terms, so t = 0 is one of the roots and what is left is a fifth degree polynomial equation with five more roots.

[Ann Kwan]

Thanks very much and yes, it would make sense to define t this way. I am wondering if there is a way to obtain an explicit solution for 5th or higher order polynomials, if such a solution exists. This will allow us to use Origin or GraphPadPrism to fit experimental data using the solution. We can use Matlab to do the fitting also but we are very new to Matlab and just starting to explore the different options.

[Roger Stafford]

@Ann: I'm afraid your chances are slim of getting a general solution to that quintic polynomial equation. If you read the article in this website:

https://en.wikipedia.org/wiki/Abel–Ruffini_theorem

you will see why. It states that "In algebra, the Abel–Ruffini theorem (also known as Abel's impossibility theorem) states that there is no general algebraic solution—that is, solution in radicals—to polynomial equations of degree five or higher with arbitrary coefficients."

[Ann Kwan]

Thanks Roger and I see. We will stick with "root" when we have 5 or more binding partners.

[Roger Stafford]

@Ann: I was mistaken in asserting that t=0 is one of the roots of the single equation in t. That is only true if 'LIDtot' is zero. Therefore this equation will in fact reduce to a sixth degree polynomial equation (with the aid of the symbolic toolbox) and will yield six roots numerically using matlab's 'roots' function.

[Ann Kwan]

@Roger: Thanks and all good. I was expecting 6 solutions (roots) for 5 competitors as there are generally n+1 solutions where n is the number of competitors. However, only one of these is meaningful in the physical experiment and we simply restrict the root value to be 0 up to LIDtot or one of the other tot (whichever appropriate, the smaller value of the two).

Answer 2

[Walter Roberson]

Z6 = (kdL4-kdM4)*(kdL3-kdL4)*(kdI2-kdL4)*(kdI1-kdL4);
Z5 = kdL4^5+(L4tot-M4tot+LIDtot-kdI1-kdI2-kdL3-kdM4-I1tot-I2tot-L3tot)*kdL4^4+((-2*kdI1-2*kdI2-2*kdL3-2*kdM4)*L4tot+(M4tot-LIDtot+kdI2+kdL3+kdM4+I2tot+L3tot)*kdI1+(M4tot-LIDtot+kdL3+kdM4+I1tot+L3tot)*kdI2+(M4tot-LIDtot+kdM4+I1tot+I2tot)*kdL3-kdM4*(LIDtot-I1tot-I2tot-L3tot))*kdL4^3+(((3*kdI2+3*kdL3+3*kdM4)*kdI1+(3*kdL3+3*kdM4)*kdI2+3*kdL3*kdM4)*L4tot+((-M4tot+LIDtot-kdL3-kdM4-L3tot)*kdI2+(-M4tot+LIDtot-kdM4-I2tot)*kdL3+kdM4*(LIDtot-I2tot-L3tot))*kdI1+((-M4tot+LIDtot-kdM4-I1tot)*kdL3+kdM4*(LIDtot-I1tot-L3tot))*kdI2+kdL3*kdM4*(LIDtot-I1tot-I2tot))*kdL4^2+((((-4*kdL3-4*kdM4)*kdI2-4*kdL3*kdM4)*kdI1-4*kdI2*kdL3*kdM4)*L4tot+(((M4tot-LIDtot+kdM4)*kdL3-kdM4*(LIDtot-L3tot))*kdI2-kdL3*kdM4*(LIDtot-I2tot))*kdI1-kdI2*kdL3*kdM4*(LIDtot-I1tot))*kdL4+5*kdM4*kdI2*(L4tot+(1/5)*LIDtot)*kdI1*kdL3;
Z4 = 6*L4tot*(((1/6)*M4tot-(1/6)*LIDtot+(1/6)*kdI1+(1/6)*kdI2+(1/6)*kdL3+(1/6)*kdM4+(1/6)*I1tot+(1/6)*I2tot+(1/6)*L3tot)*kdL4^4+(((1/6)*kdI1+(1/6)*kdI2+(1/6)*kdL3+(1/6)*kdM4)*L4tot+(-(1/3)*M4tot+(1/3)*LIDtot-(1/3)*kdI2-(1/3)*kdL3-(1/3)*kdM4-(1/3)*I2tot-(1/3)*L3tot)*kdI1+(-(1/3)*L3tot-(1/3)*M4tot+(1/3)*LIDtot-(1/3)*kdL3-(1/3)*kdM4-(1/3)*I1tot)*kdI2+(-(1/3)*I1tot-(1/3)*I2tot-(1/3)*M4tot+(1/3)*LIDtot-(1/3)*kdM4)*kdL3+(1/3)*kdM4*(LIDtot-I1tot-I2tot-L3tot))*kdL4^3+(((-(1/2)*kdI2-(1/2)*kdL3-(1/2)*kdM4)*kdI1+(-(1/2)*kdL3-(1/2)*kdM4)*kdI2-(1/2)*kdL3*kdM4)*L4tot+(((1/2)*M4tot-(1/2)*LIDtot+(1/2)*kdL3+(1/2)*kdM4+(1/2)*L3tot)*kdI2+((1/2)*M4tot-(1/2)*LIDtot+(1/2)*kdM4+(1/2)*I2tot)*kdL3-(1/2)*kdM4*(LIDtot-I2tot-L3tot))*kdI1+(((1/2)*M4tot-(1/2)*LIDtot+(1/2)*I1tot+(1/2)*kdM4)*kdL3-(1/2)*kdM4*(LIDtot-I1tot-L3tot))*kdI2-(1/2)*kdL3*kdM4*(LIDtot-I1tot-I2tot))*kdL4^2+((((kdL3+kdM4)*kdI2+kdL3*kdM4)*kdI1+kdI2*kdL3*kdM4)*L4tot+(((-(2/3)*M4tot+(2/3)*LIDtot-(2/3)*kdM4)*kdL3+(2/3)*kdM4*(LIDtot-L3tot))*kdI2+(2/3)*kdL3*kdM4*(LIDtot-I2tot))*kdI1+(2/3)*kdI2*kdL3*kdM4*(LIDtot-I1tot))*kdL4-(5/3)*kdM4*kdI2*(L4tot+(1/2)*LIDtot)*kdI1*kdL3);
Z3 = -(4*(((-(1/4)*kdM4-(1/4)*I2tot-(1/4)*L3tot-(1/4)*M4tot+(1/4)*LIDtot-(1/4)*kdI2-(1/4)*kdL3)*kdI1+(-(1/4)*kdL3-(1/4)*kdM4-(1/4)*I1tot-(1/4)*L3tot-(1/4)*M4tot+(1/4)*LIDtot)*kdI2+(-(1/4)*kdM4-(1/4)*I1tot-(1/4)*I2tot-(1/4)*M4tot+(1/4)*LIDtot)*kdL3+(1/4)*kdM4*(LIDtot-I1tot-I2tot-L3tot))*kdL4^3+(((-(1/4)*kdI2-(1/4)*kdL3-(1/4)*kdM4)*kdI1+(-(1/4)*kdL3-(1/4)*kdM4)*kdI2-(1/4)*kdL3*kdM4)*L4tot+(((3/4)*M4tot-(3/4)*LIDtot+(3/4)*kdL3+(3/4)*kdM4+(3/4)*L3tot)*kdI2+((3/4)*M4tot-(3/4)*LIDtot+(3/4)*kdM4+(3/4)*I2tot)*kdL3-(3/4)*kdM4*(LIDtot-I2tot-L3tot))*kdI1+(((3/4)*M4tot-(3/4)*LIDtot+(3/4)*I1tot+(3/4)*kdM4)*kdL3-(3/4)*kdM4*(LIDtot-I1tot-L3tot))*kdI2-(3/4)*kdL3*kdM4*(LIDtot-I1tot-I2tot))*kdL4^2+((((kdL3+kdM4)*kdI2+kdL3*kdM4)*kdI1+kdI2*kdL3*kdM4)*L4tot+(((-(3/2)*M4tot+(3/2)*LIDtot-(3/2)*kdM4)*kdL3+(3/2)*kdM4*(LIDtot-L3tot))*kdI2+(3/2)*kdL3*kdM4*(LIDtot-I2tot))*kdI1+(3/2)*kdI2*kdL3*kdM4*(LIDtot-I1tot))*kdL4-(5/2)*kdI1*kdI2*kdL3*kdM4*(L4tot+LIDtot)))*L4tot^2;
Z2 = ((((M4tot-LIDtot+kdL3+kdM4+L3tot)*kdI2+(M4tot-LIDtot+kdM4+I2tot)*kdL3-kdM4*(LIDtot-I2tot-L3tot))*kdI1+((M4tot-LIDtot+kdM4+I1tot)*kdL3-kdM4*(LIDtot-I1tot-L3tot))*kdI2-kdL3*kdM4*(LIDtot-I1tot-I2tot))*kdL4^2+((((kdL3+kdM4)*kdI2+kdL3*kdM4)*kdI1+kdI2*kdL3*kdM4)*L4tot+(((-4*M4tot+4*LIDtot-4*kdM4)*kdL3+4*kdM4*(LIDtot-L3tot))*kdI2+4*kdL3*kdM4*(LIDtot-I2tot))*kdI1+4*kdI2*kdL3*kdM4*(LIDtot-I1tot))*kdL4-5*kdI1*kdI2*kdL3*kdM4*(L4tot+2*LIDtot))*L4tot^3;
Z1 = (((((M4tot-LIDtot+kdM4)*kdL3-kdM4*(LIDtot-L3tot))*kdI2-kdL3*kdM4*(LIDtot-I2tot))*kdI1-kdI2*kdL3*kdM4*(LIDtot-I1tot))*kdL4+kdI1*kdI2*kdL3*kdM4*(L4tot+5*LIDtot))*L4tot^4;
Z0 = -L4tot^5*LIDtot*kdI1*kdI2*kdL3*kdM4;

Z = roots([Z6, Z5, Z4, Z3, Z2, Z1, Z0]);

x = -I1tot*kdL4*Z/((kdI1-kdL4)*Z-L4tot*kdI1);
y = -I2tot*kdL4*Z/((kdI2-kdL4)*Z-L4tot*kdI2);
z = -L3tot*kdL4*Z/((kdL3-kdL4)*Z-L4tot*kdL3);
w = Z;
v = (-(kdL3-kdL4)*(kdI2-kdL4)*(kdI1-kdL4)*Z^5+(-kdL4^4+(-LIDtot+kdI1+kdI2+kdL3+I1tot+I2tot+L3tot-L4tot)*kdL4^3+((2*kdI1+2*kdI2+2*kdL3)*L4tot+(LIDtot-kdI2-kdL3-I2tot-L3tot)*kdI1+(LIDtot-kdL3-I1tot-L3tot)*kdI2+kdL3*(LIDtot-I1tot-I2tot))*kdL4^2+(((-3*kdI2-3*kdL3)*kdI1-3*kdI2*kdL3)*L4tot+((-LIDtot+kdL3+L3tot)*kdI2-kdL3*(LIDtot-I2tot))*kdI1-kdI2*kdL3*(LIDtot-I1tot))*kdL4+kdI1*kdI2*kdL3*(LIDtot+4*L4tot))*Z^4-4*L4tot*(((1/4)*kdI1+(1/4)*kdI2+(1/4)*kdL3+(1/4)*I1tot-(1/4)*LIDtot+(1/4)*I2tot+(1/4)*L3tot)*kdL4^3+(((1/4)*kdI1+(1/4)*kdI2+(1/4)*kdL3)*L4tot+((1/2)*LIDtot-(1/2)*I2tot-(1/2)*L3tot-(1/2)*kdI2-(1/2)*kdL3)*kdI1+((1/2)*LIDtot-(1/2)*L3tot-(1/2)*kdL3-(1/2)*I1tot)*kdI2+(1/2)*kdL3*(LIDtot-I1tot-I2tot))*kdL4^2+(((-(3/4)*kdI2-(3/4)*kdL3)*kdI1-(3/4)*kdI2*kdL3)*L4tot+(((3/4)*kdL3-(3/4)*LIDtot+(3/4)*L3tot)*kdI2-(3/4)*kdL3*(LIDtot-I2tot))*kdI1-(3/4)*kdI2*kdL3*(LIDtot-I1tot))*kdL4+kdI1*kdI2*kdL3*(LIDtot+(3/2)*L4tot))*Z^3+(6*((((1/6)*LIDtot-(1/6)*kdI2-(1/6)*kdL3-(1/6)*I2tot-(1/6)*L3tot)*kdI1+((1/6)*LIDtot-(1/6)*kdL3-(1/6)*I1tot-(1/6)*L3tot)*kdI2+(1/6)*kdL3*(LIDtot-I1tot-I2tot))*kdL4^2+(((-(1/6)*kdI2-(1/6)*kdL3)*kdI1-(1/6)*kdI2*kdL3)*L4tot+((-(1/2)*LIDtot+(1/2)*kdL3+(1/2)*L3tot)*kdI2-(1/2)*kdL3*(LIDtot-I2tot))*kdI1-(1/2)*kdI2*kdL3*(LIDtot-I1tot))*kdL4+kdI1*kdI2*kdL3*(LIDtot+(2/3)*L4tot)))*L4tot^2*Z^2-4*L4tot^3*((((-(1/4)*LIDtot+(1/4)*kdL3+(1/4)*L3tot)*kdI2-(1/4)*kdL3*(LIDtot-I2tot))*kdI1-(1/4)*kdI2*kdL3*(LIDtot-I1tot))*kdL4+kdI1*kdI2*kdL3*(LIDtot+(1/4)*L4tot))*Z+L4tot^4*LIDtot*kdI1*kdI2*kdL3)/(((kdI1-kdL4)*Z-L4tot*kdI1)*((kdI2-kdL4)*Z-L4tot*kdI2)*((kdL3-kdL4)*Z-L4tot*kdL3)*(-L4tot+Z));

... except that you need to vectorize these, which I have not done in my transcription.

The polynomial turns out to be degree 6.

[Ann Kwan]

Thanks very much and I see. Is there a way to easily set up the equations/code so we can easily change between having different number of competitors (from 5 to 10 say, which is what we have experimentally)? For the 1 to 4 case, we can obtain the explicit solution from solve.

[Walter Roberson]

It looks to me as if it might be an eigenvalue problem??

[Ann Kwan]

Possibly. Unfortunately my maths isn't strong...