[R] R emulation of FindRoot in Mathematica

Thu Jan 19 15:42:35 CET 2023

Hi Tim  - the results from  this paper have had big impact on quite 
complicated problems of physiology so I thought I might like to 
understand how they came by, calculating them myself. Writing up the 
statements from which the inferences in the paper were generated using 
Mathematica were meant to inspire clever friends to suggest probably 
quite simple ways to handle the problem in R.

All best wishes
Troels

Den 19-01-2023 kl. 15:28 skrev Ebert,Timothy Aaron:
> This is a poster child for why we like open source software. "I dump numbers into a black box and get numbers out but I cannot verify how the numbers out were calculated so they must be correct" approach to analysis does not really work for me.
> Tim
>
> -----Original Message-----
> From: R-help <r-help-bounces using r-project.org> On Behalf Of Troels Ring
> Sent: Thursday, January 19, 2023 9:18 AM
> To: Valentin Petzel <valentin using petzel.at>; r-help mailing list <r-help using r-project.org>
> Subject: Re: [R] R emulation of FindRoot in Mathematica
>
> [External Email]
>
> Thanks,   Valentin for the suggestion. I'm not sure I can go that way. I
> include below the statements from the paper containing the knowledge on the basis of which I would like to know at specified [H] the concentration of each of the many metabolites given the constraints. I have tried to contact the author to get the full code but it seems difficult.
>
> BW Troels
>
>
> 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
>
> tmg <- mg + mgatp + mghatp + mg2atp + mgadp + mghadp + mg2adp + mgpi + kghpi + mghhpi +
>     mghpcr + mglactate
>
> tk <- k + katp + kadp + kpi + khpi + khhpi + kpcr + khpcr + khhpcr
>
> tlactate <- lactate + hlactate + mglactate
>
> # conditions
>
> tatp <- 0.008
> tpcr <- 0.042
> tcr <- 0.004
> tadp <- 0.00001
> tpi <- 0.003
> tlactate <- 0.005
>
> # free K and Mg constrained to be fixed
> #
> mg <- 0.0006
> k <- 0.12
>
> Den 19-01-2023 kl. 12:11 skrev Valentin Petzel:
>> Hello Troels,
>>
>>
>> As fair as I understand you attempt to numerically solve a system of
>> non linear equations in multiple variables in R. R does not provide
>> this functionality natively, but have you tried multiroot from the
>> rootSolve package:
>>
>>
>> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fcran
>> .r-project.org%2Fweb%2Fpackages%2FrootSolve%2FrootSolve.pdf&data=05%7C
>> 01%7Ctebert%40ufl.edu%7C7cb98cd926b34284cd5f08dafa28026c%7C0d4da0f84a3
>> 14d76ace60a62331e1b84%7C0%7C0%7C638097347110882622%7CUnknown%7CTWFpbGZ
>> sb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3
>> D%7C3000%7C%7C%7C&sdata=D9A3fwJ5x7GbEV4A01wncLUil7szTdSPul5vd0lsSBw%3D
>> &reserved=0
>>
>>
>> multiroot is called like
>>
>>
>> multiroot(f, start, ...)
>>
>>
>> where f is a function of one argument which is a vector of n values
>> (representing the n variables) and returning a vector of d values
>> (symbolising the d equations) and start is a vector of length n.
>>
>>
>> E.g. if we want so solve
>>
>>
>> x^2 + y^2 + z^2 = 1
>>
>> x^3-y^3 = 0
>>
>> x - z = 0
>>
>>
>> (which is of course equivalent to x = y = z, x^2 + y^2 + z^2 = 1, so x
>> = y = z = ±sqrt(1/3) ~ 0.577)
>>
>>
>> we'd enter
>>
>>
>> f <- function(x) c(x[1]**2 + x[2]**2 + x[3]**2 - 1, x[1]**3 - x[2]**3,
>> x[1] - x[3])
>>
>>
>> multiroot(f, c(0,0,0))
>>
>>
>> which yields
>>
>>
>> $root
>>
>> [1] 0.5773502 0.5773505 0.5773502
>>
>>
>> $f.root
>>
>> [1] 1.412261e-07 -2.197939e-07  0.000000e+00
>>
>>
>> $iter
>>
>> [1] 31
>>
>>
>> $estim.precis
>>
>> [1] 1.2034e-07
>>
>>
>> Best regards,
>>
>> Valentin
>>
>>
>> Am Donnerstag, 19. Jänner 2023, 10:41:22 CET schrieb Troels Ring:
>>
>>> Hi friends - I hope this is not a misplaced question. From the
>>> literature (Kushmerick AJP 1997;272:C1739-C1747) I have a series of
>>> Mathematica equations which are solved together to yield over
>>> different
>>> pH values the concentrations of metabolites in skeletal muscle using
>>> the
>>> Mathematica function FindRoot((E1,E2...),(V2,V2..)] where E is a
>>> list of
>>> equations and V list of variables.  Most of the equations are
>>> individual
>>> binding reactions of the form 10^6.494*atp*h == hatp and next
>>> 10^9.944*hatp*h ==hhatp describing binding of singe protons or Mg or
>>> K
>>> to ATP or creatin for example, but we also have constraints giving
>>> total
>>> concentrations of say ATP i.e. ATP + ATPH, ATPH2..ATP.Mg
>>> I have, without success, tried to find ways to do this in R - I have
>>> 36
>>> equations on 36 variables and 8 equations on total concentrations.
>>> As
>>> far as I can see from the definition of FindRoot in Wolfram, Newton
>>> search or secant search is employed.
>>> I'm on Windows R 4.2.2
>>> Best wishes
>>> Troels Ring, MD
>>> Aalborg, Denmark
>>> ______________________________________________
>>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fst
>>> at.ethz.ch%2Fmailman%2Flistinfo%2Fr-help&data=05%7C01%7Ctebert%40ufl
>>> .edu%7C7cb98cd926b34284cd5f08dafa28026c%7C0d4da0f84a314d76ace60a6233
>>> 1e1b84%7C0%7C0%7C638097347110882622%7CUnknown%7CTWFpbGZsb3d8eyJWIjoi
>>> MC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C
>>> %7C%7C&sdata=7DTBQItdQpAqK%2FCS1%2BqQvYdlvjyJMjzTOXhoS6AY%2FJQ%3D&re
>>> served=0
>>> PLEASE do read the posting guide
>> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r
>> -project.org%2Fposting-guide.html&data=05%7C01%7Ctebert%40ufl.edu%7C7c
>> b98cd926b34284cd5f08dafa28026c%7C0d4da0f84a314d76ace60a62331e1b84%7C0%
>> 7C0%7C638097347110882622%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiL
>> CJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=m0
>> 0vZP75nk9icL6H8Gc0dH1bHhkRCS9I5N27uORQmQ0%3D&reserved=0
>>
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>          [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-help&data=05%7C01%7Ctebert%40ufl.edu%7C7cb98cd926b34284cd5f08dafa28026c%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C638097347110882622%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=7DTBQItdQpAqK%2FCS1%2BqQvYdlvjyJMjzTOXhoS6AY%2FJQ%3D&reserved=0
> PLEASE do read the posting guide https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r-project.org%2Fposting-guide.html&data=05%7C01%7Ctebert%40ufl.edu%7C7cb98cd926b34284cd5f08dafa28026c%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C638097347110882622%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=m00vZP75nk9icL6H8Gc0dH1bHhkRCS9I5N27uORQmQ0%3D&reserved=0
> and provide commented, minimal, self-contained, reproducible code.