R emulation of FindRoot in Mathematica

Dear Troels,

I send you an updated version of the response. I think that a hybrid 
approach is probably the best solution:
- Input / Starting values = free: ATP, ADP, Crea, CreaP, lactate, 
inorganic phosphate;
- Output: diff(Total, given total value);
- I assume that the pH is given;
- I did NOT check all individual eqs;

library(rootSolve)

solve.AcidSpecies = function(x, H, Mg=0.0006, K = 0.12) {
??? # ... the eqs: ...;
??? ATPTotal = KaATPH * ATP * H + KaATPH2 * KaATPH * ATP * H^2 +
??? + KaATPH3 * KaATPH2 * KaATPH * ATP * H^3 + KdATPMg * ATP * Mg +
??? + KdATPHMg * ATP * H * Mg + KdATPMg2 * ATP * Mg^2 + KdATPK * ATP * K;
??? ### Output:
??? total = c(
??? ??? ATPTotal - 0.008,
??? ??? ADPTotal - 1E-5,
??? ??? CreaTotal - 0.004,
??? ??? CreaPTotal - 0.042,
??? ??? PiTotal - 0.003,
??? ??? LactateTotal - 0.005);
??? return(total);
}

KaATPH = 10^6.494; # ...
x0 = c(ATP = 0.008, ADP = 1E-5,
??? Crea = 0.004, CreaP = 0.042, Pi = 0.003, Lactate = 0.005) / 2;
x = multiroot(solve.AcidSpecies, x0, H = 4E-8);


print(x)


I think that it is possible to use the eqs directly as provided 
initially (with some minor corrections). You only need to output the 
known totals (as a diff), see full code below.


Sincerely,


Leonard



library(rootSolve)

solve.AcidSpecies = function(x, H, Mg=0.0006, k = 0.12) {
?? ?# with(list(x), { seems NOT to work with multiroot });
?? ?atp = x[1]; adp = x[2]; pi = x[3];
?? ?pcr = x[4]; cr = x[5]; lactate = x[6];

###
hatp <- 10^6.494*H*atp
hhatp <- 10^3.944*H*hatp
hhhatp <- 10^1.9*H*hhatp
hhhhatp? <- 10*H*hhhatp
mgatp <- 10^4.363*atp*Mg
mghatp <- 10^2.299*hatp*Mg
mg2atp <- 10^1-7*Mg*mgatp
katp <- 10^0.959*atp*k

hadp <- 10^6.349*adp*H
hhadp <- 10^3.819*hadp*H
hhhadp <- 10*H*hhadp
mgadp <- 10^3.294*Mg*adp
mghadp <- 10^1.61*Mg*hadp
mg2adp <- 10*Mg*mgadp
kadp <- 10^0.82*k*adp

hpi <- 10^11.616*H*pi
hhpi <- 10^6.7*H*hpi
hhhpi <- 10^1.962*H*hhpi
mgpi <- 10^3.4*Mg*pi
mghpi <- 10^1.946*Mg*hpi
mghhpi <- 10^1.19*Mg*hhpi
kpi <- 10^0.6*k*pi
khpi <- 10^1.218*k*hpi
khhpi <- 10^-0.2*k*hhpi

hpcr <- 10^14.3*H*pcr
hhpcr <- 10^4.5*H*hpcr
hhhpcr <- 10^2.7*H*hhpcr
hhhhpcr <- 100*H*hhhpcr
mghpcr <- 10^1.6*Mg*hpcr
kpcr <- 10^0.74*k*pcr
khpcr <- 10^0.31*k*hpcr
khhpcr <- 10^-0.13*k*hhpcr

hcr <- 10^14.3*H*cr
hhcr <- 10^2.512*H*hcr

hlactate <- 10^3.66*H*lactate
mglactate <- 10^0.93*Mg*lactate

tatp <- atp + hatp + hhatp + hhhatp + mgatp + mghatp + mg2atp + katp
tadp <- adp + hadp + hhadp + hhhadp + mghadp + mgadp + mg2adp + kadp
tpi? <- pi + hpi + hhpi + hhhpi + mgpi + mghpi + mghhpi + kpi + khpi + 
khhpi
tpcr <- pcr + hpcr + hhpcr + hhhpcr + hhhhpcr + mghpcr + kpcr + khpcr 
+ khhpcr
tcr? <- cr + hcr + hhcr
tlactate <- lactate + hlactate + mglactate
# tmg <- Mg + mgatp + mghatp + mg2atp + mgadp + mghadp + mg2adp + mgpi +
#??? kghpi + mghhpi + mghpcr + mglactate
# tk <- k + katp + kadp + kpi + khpi + khhpi + kpcr + khpcr + khhpcr


total = c(
?? ?tatp - 0.008,
?? ?tadp - 0.00001,
?? ?tpi - 0.003,
?? ?tpcr - 0.042,
?? ?tcr - 0.004,
?? ?tlactate - 0.005)
return(total);
# })
}

# conditions

x0 = c(
?? ?atp = 0.008,
?? ?adp = 0.00001,
?? ?pi = 0.003,
?? ?pcr = 0.042,
?? ?cr = 0.004,
?? ?lactate = 0.005
) / 3;
# tricky to get a positive value !!!
x0[1] = 0.001; # still NOT positive;

x = multiroot(solve.AcidSpecies, x0, H = 4E-8)


On 1/23/2023 12:37 AM, Leonard Mada wrote:

Dear Troels,

The system that you mentioned needs to be transformed first. The 
equations are standard acid-base equilibria-type equations in 
analytic chemistry.

ATP + H <-> ATPH
ATPH + H <-> ATPH2
ATPH2 + H <-> ATPH3
[...]
The total amount of [ATP] is provided, while the concentration of the 
intermediates are unknown.

Q.) It was unclear from your description:
Do you know the pH?
Or is the pH also unknown?

I believe that the system is exactly solvable. The "multivariable" 
system/solution may be easier to write down: but is uglier to solve, 
as the "system" is under-determined. You can use optim in such cases, 
see eg. an example were I use it:
https://github.com/discoleo/R/blob/master/Stat/Polygons.Examples.R


a2 = optim(c(0.9, 0.5), polygonOptim, d=x)
# where the function polygonOptim() computes the distance between the 
starting-point & ending point of the polygon;
# (the polygon is defined only by the side lengths and optim() tries 
to compute the angles);
# optimal distance = 0, when the polygon is closed;
# Note: it is possible to use more than 2 starting values in the 
example above (the version with optim) works quit well;
# - but you need to "design" the function that is optimized for your 
particular system, e.g.
#?? by returning: c((ATPTotal - value)^2, (ADPTotal - value)^2, ...);


S.1.) Exact Solution
ATP system: You can express all components as eqs of free ATP, [ATP], 
and [H], [Mg], [K].
ATPH = KaATPH * ATP * H;
ATPH2 = KaATPH2 * ATPH * H
= KaATPH2 * KaATPH * ATP * H^2;
[...]

Then you substitute these into the total-eqs. It is a bit uglier, as 
you have 5 coupled systems: ATP, ADP, Creatinine, CreaP and Lactate. 
This is why I thought that it may be easier to do it, at least 
partially/hybrid, in a "multivariate" way.

# Total ATP:
ATPTotal = KaATPH * ATP * H + KaATPH2 * KaATPH * ATP * H^2 +
??? + KaATPH3 * KaATPH2 * KaATPH * ATP * H^3 + ...;

Note:
- the system above is solvable, if the pH is known; it may be 
undetermined if the pH is NOT known (but I am not very certain: 
easier to spot only when all the eqs are written down nicely);
- the eqs for total K & total Mg: are not useful to solve the system, 
as there are NO constraints on the total values; they will yield only 
some numerical values (which may be interesting for the analysis);
- total K & total Mg can be computed after solving the system;

Please let me know if you need further assistance.

Sincerely,

Leonard