[R] R emulation of FindRoot in Mathematica

[Troels Ring]

Dear Leonard - thanks a lot for the solution! Yes, you are right: pH is 
presumed known, so the effort is repeated at a series of chosen pH 
values.  It appears that no matter x0 the estimate for atp is always 
-6.239560e-01 which even accounting for offset of 0.008 is negative 
which is not possible and the estimate of the other species also appear 
to be constant no matter x0? The algorithm needs only 2 or 3 iterations 
and doesn't complain. I have been through the equations and they appear 
correct. Wonder what I did wrong?

Below some output

Best wishes
Troels

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)
x
# $root
# atp           adp            pi           pcr cr       lactate
# -6.239560e-01  3.254998e-06  5.581774e-08  4.142785e-09 5.011807e-10  
4.973691e-03

x0 <- 1000*x0
x = multiroot(solve.AcidSpecies, x0, H = 4E-8)
x
# $root
# atp           adp            pi           pcr cr       lactate
# -6.239560e-01  3.254998e-06  5.581774e-08  4.142785e-09 5.011807e-10  
4.973691e-03

x0 <- x0*1e-6
x = multiroot(solve.AcidSpecies, x0, H = 4E-8)
x
# $root
# atp           adp            pi           pcr cr       lactate
# -6.239560e-01  3.254998e-06  5.581774e-08  4.142785e-09 5.011807e-10  
4.973691e-03


Den 23-01-2023 kl. 01:24 skrev Leonard Mada:
> 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
>>
>>
	[[alternative HTML version deleted]]