[R] R emulation of FindRoot in Mathematica

[Leonard Mada]

Dear Troels,


There might be an error in one of the eqs:

# [modified] TODO: check;
mg2atp <- 10^(-7)*Mg*mgatp;

This version works:
x0 = c(
     atp = 0.008,
     adp = 0.00001,
     pi = 0.003,
     pcr = 0.042,
     cr = 0.004,
     lactate = 0.005
) / 30;
# solved the positive value
x0[1] = 1E-6;

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

# Results:
#         atp          adp           pi          pcr cr      lactate
# 4.977576e-04 3.254998e-06 5.581774e-08 4.142785e-09 5.011807e-10 
4.973691e-03


Sincerely,


Leonard


On 1/23/2023 2:24 AM, Leonard Mada wrote:
> 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]]