[CHNOSZ-commits] r940 - in pkg/CHNOSZ: . R inst inst/tinytest man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Dec 7 05:49:46 CET 2025
Author: jedick
Date: 2025-12-07 05:49:45 +0100 (Sun, 07 Dec 2025)
New Revision: 940
Added:
pkg/CHNOSZ/R/phosphorylate.R
pkg/CHNOSZ/inst/tinytest/test-phosphorylate.R
pkg/CHNOSZ/man/phosphorylate.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/inst/NEWS.Rd
Log:
Calculate affinity of phosphorylation reactions
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2025-12-06 09:59:54 UTC (rev 939)
+++ pkg/CHNOSZ/DESCRIPTION 2025-12-07 04:49:45 UTC (rev 940)
@@ -1,6 +1,6 @@
-Date: 2025-12-06
+Date: 2025-12-07
Package: CHNOSZ
-Version: 2.2.0-6
+Version: 2.2.0-7
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2025-12-06 09:59:54 UTC (rev 939)
+++ pkg/CHNOSZ/NAMESPACE 2025-12-07 04:49:45 UTC (rev 940)
@@ -60,7 +60,9 @@
# added 20220416
"rank.affinity",
# added 20220620
- "stack_mosaic"
+ "stack_mosaic",
+# added 20251207
+ "phosphorylate"
)
# Load shared objects
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2025-12-06 09:59:54 UTC (rev 939)
+++ pkg/CHNOSZ/R/basis.R 2025-12-07 04:49:45 UTC (rev 940)
@@ -96,7 +96,10 @@
# We don't accept any of those
if(is.list(ispecies)) ina <- ina | sapply(ispecies,length) > 1
}
- if(any(ina)) stop(paste("species not available:", paste(species[ina], "(", state[ina], ")", sep = "", collapse = " ")))
+ if(any(ina)) {
+ if(is.null(state)) stop(paste("species not available:", paste(species[ina], sep = "", collapse = " ")))
+ else stop(paste("species not available:", paste(species[ina], "(", state[ina], ")", sep = "", collapse = " ")))
+ }
# Are we adding a species to an existing basis definition? 20220208
if(add) {
Added: pkg/CHNOSZ/R/phosphorylate.R
===================================================================
--- pkg/CHNOSZ/R/phosphorylate.R (rev 0)
+++ pkg/CHNOSZ/R/phosphorylate.R 2025-12-07 04:49:45 UTC (rev 940)
@@ -0,0 +1,190 @@
+# Calculate affinity of phosphorylation reactions taking account of speciation
+# 20251206 first version (extracted from sugars paper script) jmd
+phosphorylate <- function(reactant, P_source, loga_reactant = 0, loga_product = 0, loga_P_source = 0, loga_P_remainder = 0, const_pH = 7, ...) {
+
+ # Affinity is calculated by summing:
+ # 1) reactant + P = product + H2O (m1)
+ # 2) P_source + H2O = P_remainder + P (m2)
+ # NOTE: reaction 2 is defined in reverse (to form H2O), then the opposite of the affinity is used for the sum
+
+ # Terminology:
+ # Complex species are species that take multiple forms (e.g. ionization or oxidation states) depending on pH or other variables
+ # A basic reaction has complex species that are valid basis species (independent components)
+ # --> This is what mosaic() deals with
+ # An extended reaction has complex species all of which cannot be used as basis species (they are non-independent)
+ # However, the energetics of an extended reaction can be modeled as a sum of basic reactions
+ # --> This script has a function for this called extended_mosaic()
+
+ # Examples:
+ # Basic reaction: acetic acid + P = acetylphosphate + H2O
+ # Extended reaction: acetic acid + PP = acetylphosphate + P
+ # Complex species: acetic acid = [acetic acid + acetate]
+
+ ## Setup reaction 1
+ if(reactant == "acetic acid") {
+ # Basic reaction: acetic acid + P = acetylphosphate + H2O
+ # Load initial species for mosaic reaction (uncharged species)
+ basis(c("acetic acid", "H3PO4", "acetylphosphate0", "O2", "H+"))
+ # The basis species we will swap through for mosaic
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("acetylphosphate0", "acetylphosphate-1", "acetylphosphate-2", "acetylphosphate-3"),
+ c("acetic acid", "acetate")
+ )
+ } else if(reactant == "glycerol") {
+ # Basic reaction: glycerol + P = 1-glycerolphosphate + H2O
+ basis(c("glycerol", "H3PO4", "3-phosphoglycerol", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("3-phosphoglycerol", "3-phosphoglycerol-1", "3-phosphoglycerol-2")
+ )
+ } else if(reactant == "adenosine_to_AMP") {
+ # Basic reaction: adenosine + P = AMP + H2O
+ basis(c("adenosine", "H3PO4", "H2AMP", "N2", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("H2AMP", "HAMP-", "AMP-2")
+ )
+ } else if(reactant == "adenosine_to_cAMP") {
+ # Basic reaction: adenosine + P = cAMP + H2O
+ basis(c("adenosine", "H3PO4", "cyclic-HAMP", "N2", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("cyclic-HAMP", "cyclic-AMP-1")
+ )
+ } else if(reactant == "ribose") {
+ # Basic reaction: ribose + P = ribose-5-phosphate + H2O
+ basis(c("ribose", "H3PO4", "ribose-5-phosphate", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("ribose-5-phosphate", "ribose-5-phosphate-1", "ribose-5-phosphate-2")
+ )
+ } else if(reactant == "uridine") {
+ # Basic reaction: uridine + P = UMP + H2O
+ basis(c("uridine", "H3PO4", "H2UMP", "N2", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("H2UMP", "HUMP-", "UMP-2")
+ )
+ } else if(reactant == "AMP") {
+ # Basic reaction: AMP + P = ADP + H2O
+ basis(c("H2AMP", "H3PO4", "H3ADP", "N2", "O2", "H+"))
+ bases <- list(
+ c("H2AMP", "HAMP-", "AMP-2"),
+ c("H3ADP", "H2ADP-", "HADP-2", "ADP-3")
+ )
+ } else if(reactant == "ADP") {
+ # Basic reaction: ADP + P = ATP + H2O
+ basis(c("H3ADP", "H3PO4", "H4ATP", "N2", "O2", "H+"))
+ bases <- list(
+ c("H3ADP", "H2ADP-", "HADP-2", "ADP-3"),
+ c("H4ATP", "H3ATP-", "H2ATP-2", "HATP-3", "ATP-4")
+ )
+ } else if(reactant == "glucose") {
+ # Basic reaction: glucose + P = glucose-6-phosphate + H2O
+ basis(c("glucose", "H3PO4", "glucose-6-phosphate", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("glucose-6-phosphate", "glucose-6-phosphate-1", "glucose-6-phosphate-2")
+ )
+ } else if(reactant == "pyruvic acid") {
+ # Basic reaction: pyruvic acid + P = phosphoenolpyruvate + H2O
+ basis(c("pyruvic acid", "H3PO4", "phosphoenolpyruvate", "O2", "H+"))
+ bases <- list(
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3"),
+ c("phosphoenolpyruvate", "phosphoenolpyruvate-1", "phosphoenolpyruvate-2", "phosphoenolpyruvate-3"),
+ c("pyruvic acid", "pyruvate")
+ )
+ } else {
+ stop(paste("unrecognized reactant:", reactant))
+ }
+
+ # Set activity of H3PO4
+ # Basic reaction: H3PO4 is the P_source
+ if(P_source == "P") basis("H3PO4", loga_P_source)
+ # Extended reaction with PP: H3PO4 is the P_remainder
+ if(P_source == "PP") basis("H3PO4", loga_P_remainder)
+ # Other extended reactions [don't set the H3PO4 activity]
+ # There is no H3PO4 in the overall reaction, but is present in the half reactions.
+ # For correct cancellation, we just use the default activity of H3PO4.
+
+ # NOTE: The initial basis definition uses the name, but setting the activity uses the elemental formula
+ formula_product <- rownames(basis())[3]
+ basis(formula_product, loga_product)
+ basis("pH", const_pH)
+ # This is used for changing the activity of the reactant
+ formula_reactant <- rownames(basis())[1]
+ basis(formula_reactant, loga_reactant)
+
+ # Generate overall reaction
+ # Don't use species("acetylphosphate0"), because that reaction is just acetylphosphate0 = acetylphosphate0
+ # Instead, form H2O from the basis species to generate the overall reaction
+ species("H2O")
+
+ # Calculate affinity of forming product from predominant basis species
+ m1 <- mosaic(bases, ...)
+ # For cAMP, double the reaction to consume one H3PO4 20251201
+ if(reactant == "adenosine_to_cAMP") m1$A.species$values <- lapply(m1$A.species$values, "*", 2)
+
+ ## Setup reaction 2
+ if(P_source == "P") {
+ # We just want the first reaction (the basic reaction)
+ m2 <- NULL
+ a12 <- m1$A.species$values[[1]]
+ P_reaction <- NULL
+ } else {
+ # Save the existing basis and species definition (used for m1)
+ basis_1 <- basis()
+ species_1 <- species()
+ # Mosaic for opposite of reaction 2:
+ if(P_source == "PP") {
+ # Form H2O from H4P2O7 and H3PO4
+ basis(c("H4P2O7", "H3PO4", "O2", "H+"))
+ basis("H4P2O7", loga_P_source)
+ basis("H3PO4", loga_P_remainder)
+ bases <- list(
+ c("H4P2O7", "H3P2O7-", "H2P2O7-2", "HP2O7-3", "P2O7-4"),
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3")
+ )
+ } else if(P_source == "acetylphosphate") {
+ # Form H2O from acetylphosphate and acetic acid
+ basis(c("acetylphosphate0", "acetic acid", "H3PO4", "O2", "H+"))
+ basis("C2H5O5P", loga_P_source)
+ basis("C2H4O2", loga_P_remainder)
+ bases <- list(
+ c("acetylphosphate0", "acetylphosphate-1", "acetylphosphate-2", "acetylphosphate-3"),
+ c("acetic acid", "acetate"),
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3")
+ )
+ } else if(P_source == "ATP") {
+ # Form H2O from ADP and ATP
+ basis(c("H4ATP", "H3ADP", "H3PO4", "N2", "O2", "H+"))
+ basis("C10H16N5O13P3", loga_P_source)
+ basis("C10H15N5O10P2", loga_P_remainder)
+ bases <- list(
+ c("H4ATP", "H3ATP-", "H2ATP-2", "HATP-3", "ATP-4"),
+ c("H3ADP", "H2ADP-", "HADP-2", "ADP-3"),
+ c("H3PO4", "H2PO4-", "HPO4-2", "PO4-3")
+ )
+ }
+ species("H2O")
+ basis("pH", const_pH)
+ # Use the bases defined above and take the other arguments (e.g. pH, T, and P) from m1
+ m_args <- m1$args
+ m_args[["bases"]] <- bases
+ # TODO: handle ionic strength
+ #if(IS > 0) m_args[["IS"]] <- IS
+ m2 <- do.call(mosaic, m_args)
+ # The affinity of the extended reaction
+ a12 <- m1$A.species$values[[1]] - m2$A.species$values[[1]]
+ # Return the formation reaction to use in describe.reaction() 20251201
+ P_reaction <- species()
+ # Restore the previous basis and species definitions for later calculations
+ thermo(basis = basis_1)
+ thermo(species = species_1)
+ }
+
+ # Put together the output
+ list(m1 = m1, m2 = m2, a12 = a12, P_reaction = P_reaction)
+
+}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2025-12-06 09:59:54 UTC (rev 939)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2025-12-07 04:49:45 UTC (rev 940)
@@ -15,7 +15,7 @@
\newcommand{\Cp}{\ifelse{latex}{\eqn{C_P}}{\ifelse{html}{\out{<I>C<sub>P</sub></I>}}{Cp}}}
\newcommand{\DG0}{\ifelse{latex}{\eqn{{\Delta}G^{\circ}}}{\ifelse{html}{\out{Δ<I>G</I>°}}{ΔG°}}}
-\section{Changes in CHNOSZ version 2.2.0-4 (2025-08-21)}{
+\section{Changes in CHNOSZ version 2.2.0-7 (2025-12-07)}{
\itemize{
@@ -31,6 +31,9 @@
\item Add stibnite (Sb\s{2}S\s{3}) from
\href{https://doi.org/10.3133/b2131}{Robie and Hemingway (1995)}.
+ \item Add \code{phosphorylate()} to calculate affinity of phosphorylation
+ reactions with pH-dependent speciation.
+
}
}
Added: pkg/CHNOSZ/inst/tinytest/test-phosphorylate.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-phosphorylate.R (rev 0)
+++ pkg/CHNOSZ/inst/tinytest/test-phosphorylate.R 2025-12-07 04:49:45 UTC (rev 940)
@@ -0,0 +1,35 @@
+# 20251207
+
+# First add thermodynamic data for sugar phosphates (from Table 3.2 of Alberty, 2003)
+mod.OBIGT("glucose-6-phosphate-2", formula = "C6H11O9P-2", G = -1763940)
+mod.OBIGT("glucose-6-phosphate-1", formula = "C6H12O9P-", G = -1800590)
+# Alberty (2003) doesn't have ΔG° for neutral glucose-6-phosphate,
+# so we calculate it from pKa1 = 1.5 (Degani and Halmann, 1966)
+DG0_G6P <- -1800590 + convert(1.5, "G")
+mod.OBIGT("glucose-6-phosphate", formula = "C6H13O9P", G = DG0_G6P)
+
+# Calculate affinity at pH 5-9 with unit activities (loga = 0)
+pH <- 5:9
+result <- phosphorylate("glucose", "ATP", pH = pH)
+
+# Extract the overall affinity
+A <- result$a12
+
+# Convert affinity to Gibbs energy in kJ/mol
+TK <- convert(25, "K")
+deltaG_J <- convert(A, "G", T = TK)
+deltaG_kJ <- deltaG_J / 1000
+
+info <- "ΔG°' becomes more negative at higher pH"
+expect_true(all(diff(deltaG_kJ) < 0), info = info)
+
+info <- "ΔG°' has expected value at pH 7"
+expect_equal(round(deltaG_kJ[3], 2), -32.15, info = info)
+
+info <- "Using const_pH argument gives expected value"
+result <- phosphorylate("glucose", "ATP", const_pH = 7)
+A <- result$a12
+TK <- convert(25, "K")
+deltaG_J <- convert(A, "G", T = TK)
+deltaG_kJ <- deltaG_J / 1000
+expect_equal(round(as.numeric(deltaG_kJ), 2), -32.15, info = info)
Added: pkg/CHNOSZ/man/phosphorylate.Rd
===================================================================
--- pkg/CHNOSZ/man/phosphorylate.Rd (rev 0)
+++ pkg/CHNOSZ/man/phosphorylate.Rd 2025-12-07 04:49:45 UTC (rev 940)
@@ -0,0 +1,121 @@
+\name{phosphorylate}
+\alias{phosphorylate}
+\title{Calculate affinity of phosphorylation reactions with pH-dependent speciation}
+\description{
+ Calculate the affinity and Gibbs energy of phosphorylation reactions accounting for pH-dependent speciation of reactants, products, and phosphate sources using the mosaic approach.
+}
+
+\usage{
+phosphorylate(reactant, P_source, loga_reactant = 0, loga_product = 0,
+ loga_P_source = 0, loga_P_remainder = 0, const_pH = 7, ...)
+}
+\arguments{
+ \item{reactant}{
+ character, the compound to be phosphorylated.
+ Options include \code{"acetic acid"}, \code{"glycerol"}, \code{"adenosine_to_AMP"}, \code{"adenosine_to_cAMP"},
+ \code{"ribose"}, \code{"uridine"}, \code{"AMP"}, \code{"ADP"}, \code{"glucose"}, or \code{"pyruvic acid"}.}
+ \item{P_source}{
+ character, the source of phosphate.
+ Use \code{"P"} for basic reactions with inorganic phosphate,
+ or \code{"PP"}, \code{"acetylphosphate"}, or \code{"ATP"} for extended reactions with other phosphate donors.}
+ \item{loga_reactant}{numeric, logarithm of activity of the reactant. Default: 0}
+ \item{loga_product}{numeric, logarithm of activity of the product. Default: 0}
+ \item{loga_P_source}{numeric, logarithm of activity of the phosphate source. Default: 0}
+ \item{loga_P_remainder}{numeric, logarithm of activity of the phosphate remainder (e.g., Pi when PP is the source). Default: 0}
+ \item{const_pH}{numeric, pH value for the calculation. Default: 7}
+ \item{\dots}{
+ additional arguments passed to \code{\link{mosaic}},
+ such as \code{pH}, \code{T} (temperature), or \code{P} (pressure) ranges for multi-dimensional calculations.}
+}
+
+\details{
+ This function calculates the affinity of phosphorylation reactions by summing the affinities of component reactions
+ while accounting for pH-dependent speciation of all complex species (reactants, products, and phosphate forms).
+ Complex species are those that exist in multiple ionization states depending on pH.
+ The \code{\link{mosaic}} function is used to calculate affinities accounting for the predominant species at each pH condition.
+
+ The calculation method depends on whether the reaction is \dQuote{basic} or \dQuote{extended}:
+
+ \strong{Basic reactions} (\code{P_source = "P"}): A single reaction of the form:
+ \preformatted{ reactant + Pi = product + H2O}
+
+ \strong{Extended reactions} (other \code{P_source} values): The sum of two reactions:
+ \preformatted{ (1) reactant + Pi = product + H2O
+ (2) P_source + H2O = P_remainder + Pi}
+
+ The overall affinity is calculated as A_total = A1 - A2, where the negative sign accounts for the reverse direction of reaction 2.
+
+ Use \code{const_pH} to perform calculations at constant pH (when \code{pH} is not provided as a range in \code{...}).
+
+ If pressure (\code{P}) is one of the arguments in \code{...},
+ then call this function with \code{P_source} as a named argument to avoid conflicting assignments due to partial argument matching.
+
+ The function returns affinity values that can be converted to standard transformed Gibbs energy using \code{\link{convert}}.
+}
+\section{Warning}{
+ This function depends on thermodynamic data for phosphorylated species that have not yet been added to OBIGT (2025-12-07).
+ The data for glucose-6-phosphate used in the examples are provisional and are known to produce inaccurate results for the Gibbs energy of phosphorylation
+ when used with the data in OBIGT for glucose, ATP, and ADP.
+}
+
+\value{
+ A list with the following components:
+ \item{m1}{the \code{\link{mosaic}} output for reaction 1 (reactant + Pi = product + H2O)}
+ \item{m2}{the \code{\link{mosaic}} output for reaction 2 (NULL for basic reactions with \code{P_source = "P"})}
+ \item{a12}{
+ numeric, the overall affinity (dimensionless, A/2.303RT).
+ For basic reactions this equals the affinity from m1; for extended reactions this is the sum of affinities from m1 and m2.}
+ \item{P_reaction}{the formation reaction definition for the P_source (NULL for basic reactions)}
+}
+\seealso{
+ \code{\link{mosaic}}, \code{\link{affinity}}, \code{\link{convert}}
+}
+
+\examples{
+\dontshow{data(thermo)}
+# Calculate standard transformed Gibbs energy (DeltaG°') for
+# glucose + ATP = glucose-6-phosphate + ADP at pH 7, 25 °C
+
+# First add thermodynamic data for sugar phosphates (from Table 3.2 of Alberty, 2003)
+mod.OBIGT("glucose-6-phosphate-2", formula = "C6H11O9P-2", G = -1763940)
+mod.OBIGT("glucose-6-phosphate-1", formula = "C6H12O9P-", G = -1800590)
+# Alberty (2003) doesn't have DeltaG° for neutral glucose-6-phosphate,
+# so we calculate it from pKa1 = 1.5 (Degani and Halmann, 1966)
+DG0_G6P <- -1800590 + convert(1.5, "G")
+mod.OBIGT("glucose-6-phosphate", formula = "C6H13O9P", G = DG0_G6P)
+
+# Calculate affinity at pH 7 with unit activities (loga = 0)
+result <- phosphorylate("glucose", "ATP", const_pH = 7)
+
+# Extract the overall affinity
+A <- result$a12
+
+# Convert affinity to Gibbs energy in kJ/mol
+TK <- convert(25, "K")
+deltaG_J <- convert(A, "G", T = TK)
+deltaG_kJ <- deltaG_J / 1000
+
+print(paste("DeltaG°' at pH 7:", round(deltaG_kJ, 2), "kJ/mol"))
+
+# Example with non-standard activities
+# log(a) = -3 for glucose, ATP, and products
+result2 <- phosphorylate("glucose", "ATP",
+ loga_reactant = -5,
+ loga_product = -4,
+ loga_P_source = -3,
+ loga_P_remainder = -2,
+ const_pH = 7)
+A2 <- result2$a12
+deltaG2_kJ <- convert(A2, "G", T = convert(25, "K")) / 1000
+print(paste("DeltaG at pH 7 with non-unit activities:",
+ round(deltaG2_kJ, 2), "kJ/mol"))
+}
+
+\references{
+ Alberty RA (2003) \emph{Thermodynamics of Biochemical Reactions}. John Wiley & Sons, Hoboken, New Jersey, 397 p. \url{https://www.worldcat.org/oclc/51242181}
+
+ Degani C and Halmann M. (1966) Solvolysis of phosphoric acid esters. Hydrolysis of glucose 6-phosphate. Kinetic and tracer studies.
+ \emph{J Am Chem Soc} \bold{88}:4075-4082. \doi{10.1021/ja00969a033}
+}
+
+\concept{Thermodynamic calculations}
More information about the CHNOSZ-commits
mailing list