[R] R emulation of FindRoot in Mathematica

Thu Jan 19 15:28:53 CET 2023

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.