I'm trying to solve an equilibrium problem for the following competing chemical reaction:
AP + P >< APa@P (K1 = APa@P/(AP*P) = 1e5)
AP + P >< APp@P (K2 = APp@P/(AP*P) = 1e4)
That is; AP can bind in two ways to P and does so with two different association constants (K1 and K2). I will therefore have 4 species in solution: unbound AP, unbound P, APa@P and APp@P. Knowing the starting concentrations AP0 and P0, fsolve should easily be able to solve all the concentrations by solving the following four nonlinear equations:
P0  P  APa@P  APp@P = 0
AP0  APa@P  APp@P = 0
APa@P/(P*AP)  K1 = 0
APp@P/(P*AP)  K2 = 0
My function file becomes:
function [Concentrations] = svdConcfunc(D)
global AP0 P0 K1 K2
Concentrations = [P0  D(2)  D(3)  D(4);
AP0  D(1)  D(3)  D(4);
D(3)./(D(1).*D(2))  K1;
D(4)./(D(1).*D(2))  K2];
return
My script, calling the above function looks like this:
clear all
close all
format long
global conc AP0 P0 K1 K2
options = optimset ('Tolfun', 1e30, 'TolX', 1e20, 'MaxFunEvals', 1e20, 'MaxIter', 1e20);
conc=load('AP_titration_conc_vector.txt')./(1e6);
P0 = 2e6;
K1 = 1e5;
K2 = 1e4;
for i = 1:length(conc);
AP0 = conc(i);
%format = [AP P APa@P APp@P ];
Guessconc = [AP0 P0 AP0 AP0 ];
CalcConc(i,:) = fsolve(@svdConcfunc, Guessconc, options);
test_should_be(i) = P0 + AP0;
end
test = sum(CalcConc, 2)
test_should_be = test_should_be'
the loaded "conc" becomes:
conc=
1.0e004 *
0
0.000800000000000
0.001600000000000
0.006380000000000
0.014300000000000
0.030000000000000
0.053000000000000
0.082800000000000
0.126000000000000
0.180000000000000
0.244000000000000
0.315000000000000
0.396000000000000
0.493000000000000
The output I get is the following:
test =
1.0e004 *
0.020000000000000
0.020656586717874
0.021314855866978
0.025281807267153
0.031972557259257
0.045604844303509
0.066255477286183
0.093836312858125
0.134837016836453
0.187048378261830
0.249671527576348
0.319655268142224
0.399862134556760
0.496206019001146
test_should_be =
1.0e004 *
0.020000000000000
0.020800000000000
0.021600000000000
0.026380000000000
0.034300000000000
0.050000000000000
0.073000000000000
0.102800000000000
0.146000000000000
0.200000000000000
0.264000000000000
0.335000000000000
0.416000000000000
0.513000000000000
So the sum of all calculated concentrations in solution is not the same (there is a difference between "test_should_be" and "test") as the sum of what I put in, which is physically impossible and I cannot for the life of me understand why it won't work. I have been able to solve far more complicated equilibrium problems (more equations, species and binding constants) before so this should be a walk in the park.
I have tried to play around with starting guesses and tolerance values but I can't seem to get a better solution.
Any help on this would be much appreciated!
/Jesper
