Energy-based bond graph models of glucose transport with SLC transporters

1Introduction¶

We presented the bond graph BG models for the solute carrier (SLC) family members SLC2A2 and SLC5A1 in Hunter et al. (2025). The models have been implemented using CellML Cuellar et al., 2003 and the model implementation is available in the Physiome Model Repository (PMR) Yu et al., 2011 at: https://models.physiomeproject.org/workspace/b65. In that workspace (and in the accompanying OMEX archive) the folders Facilitated transporter and Electrogenic cotransporter hold the models of SLC2A2 and SLC5A1, respectively. Brief descriptions of the CellML model files can be found in the PMR exposure: https://models.physiomeproject.org/e/cd3/. In Sections 3 & 4 we provide detailed introductions to these two exemplar models from Hunter et al. (2025). The instructions for reproducing all the simulation experiments presented in Hunter et al. (2025) are provided in Sections 3.3 & 4.2 for SLC2A2 and SLC5A1, respectively, with instructions on how to initialize these simulation experiments given in Section 2.

2Model and simulation setup¶

In addition to the CellML models, this study makes use of Python scripts to facilitate simulation and presentation. All files required to reproduce and reuse this study can be obtained by downloading the associated OMEX archive. Alternatively, readers familiar with git may clone the repository at https://models.physiomeproject.org/workspace/b65 using their preferred git client. Please note that as the repository uses git submodules, cloning should be performed recursively. With Git version 2.13 and later, use git clone --recurse-submodules or follow the appropriate instructions based on your Git client.

After obtaining the required files, a Python environment with the required capabilities can be created by installing the dependencies listed in the requirements.txt file. This can be done using a command similar to pip install -r requirements.txt or by following the appropriate procedure for the chosen platform.

The src folder contains Python scripts for running simulations, processing simulation results, and plotting data. The primary scripts relevant to this manuscript are summarised here as well as mentioned as they are used in the following sections.

• sim_GLUT2.py: Performs simulations of SLC2A2 and saves the results to </Facilitated transporter/CellMLV2/sim_results>.

• sim_SGLT1.py: Performs simulations of SLC5A1 and saves the results to </Electrogenic cotransporter/CellMLV2/sim_results>.

• mergeData_GLUT2.py: Prepares the data in </Facilitated transporter/CellMLV2/sim_results> for plotting.

• mergeData_SGLT1.py: Prepares the data in </Electrogenic cotransporter/CellMLV2/sim_results> for plotting.

• plot_GLUT2.py: Plots the simulation results for SLC2A2.

• plot_SGLT1.py: Plots the simulation results for SLC5A1.

Other Python scripts were used to generate the SED-ML files in </Facilitated transporter/CellMLV2/> and </Electrogenic cotransporter/CellMLV2/>. These SED-ML files are provided so users do not need to run the scripts again.

3SLC2A2 bond graph model parameterisation¶

The SLC2A2 (protein name GLUT2) uses the extracellular to intracellular glucose concentration gradient to drive transmembrane transport of glucose in a process called ‘facilitated diffusion’, and we replicated the bond graph diagram in Figure 1 for convenience. The kinetic data that we used to obtain the parameters of the bond graph model were from Lowe & Walmsley (1986), and the kinetic model diagram is shown in Figure 2. Note that the notation and the parameter names in the kinetic diagram are different from the bond graph. Additionally, the bond graph model of SLC2A2 can be generalised and parameterised to represent any member of the SLC2 family. The data from Lowe & Walmsley (1986) were measured for SLC2A1/GLUT1, while SLC2A2/GLUT2 is used in this Physiome paper to remain consistent with the primary reference Hunter et al. (2025).

Table 1 lists the kinetic parameters, which are the rate constants associated with the forward and reverse reaction fluxes in the traditional mass-action equations. The first column of Table 1 is the corresponding parameter names for the reactions in the bond graph where the subscript indicates the reaction number, while the second column lists the kinetic parameter names in Lowe & Walmsley (1986).

3.1Bond graph parameters¶

Pan (2019) introduced a method to convert kinetic parameters to bond graph parameters provided that the kinetic models are thermodynamically consistent. Here, we use the SLC2A2 example to detail the link between the kinetic parameters in traditional mass-action equations and bond graph parameters. In particular, kinetic models often capture the fluxes $J$ using the unit such as $mM.s^{-1}$, while bond graph usually explicitly models the flow rate $v$ using the unit such as $fmol.s^{-1}$ and potentials $\mu$ using the unit $J.mol^{-1}$ for biochemical reactions. When we convert the kinetic parameters to bond graph parameters, we need to consider such dimensional differences.

Conventionally, the rate of biochemical reactions can be described by the law of mass action. The rate of the forward reaction $J^+$ ($mM.s^{-1}$) (also known as forward flux) is proportional to the amount of the reactants (E, as shown in Equation 1), while the rate of the reverse reaction $J^-$ ($mM.s^{-1}$) (also known as reverse flux) is proportional to the amount of the products (P, as shown in Equation 2). $[C_i]$ ($mM$) and $[C_j]$ ($mM$) are the concentrations of reactants and products, $k^+$ and $k^-$ are the forward and reverse rate constants, while $\nu_i^f$ and $\nu_j^r$ are the corresponding stoichiometric coefficients.

The total rate of reaction $J~$ ($mM.s^{-1}$) (i.e., net flux) is expressed in Equation 3.

In the case of the reaction network of SLC2A2 (Figure 1), we can express the fluxes using Equation 4.

The bond graph formulation highlights thermodynamic consistency and the flow of chemical species is driven by the chemical potentials Pan, 2019. The chemical potential $\mu_i$ of a specicies $i$ is determined by the molar amount of $q_i$ ($fmol$), shown in Equation 5, where $R=8.314$ ($J.K^{-1}.mol^{-1}$) is the ideal gas constant, $T$ ($K$) is the absolute temperature and $K_i$ ($fmol^{-1}$) is the thermodynamic constant of the species.

The rate of a reaction $v_R$ ($fmol.s^{-1}$) can be expressed using the Marcelin-de Donder equation (Equation 6). $A_R^f$ ($J.mol^{-1}$) is the forward affinity (the total chemical potential of the reactants) and $A_R^r$ ($J.mol^{-1}$) is the reverse affinity (the total chemical potential of the products).

For example, the flow rate of the reaction $Re_3$ in Figure 1 can be given by Equations 7, 8 and 9.

Equation 7 can be rearranged as Equation 10 by substituting Equations 8 and 9 into Equations 7.

The flow rates of the reaction network of SLC2A2 (Figure 1), can be expressed using Equation 11.

The relationship between the molar amount $q_i$ ($fmol$) and the concentration $[C_i]$ ($mM$) of a species is $q_i=[C_i]V_i$, where $V_i$ ($pL$) is the volume of the compartment in which the species resides. Equation 11 can be rewritten as Equation 12 if we incorporate the concentrations of species rather than the molar amount.

By comparing Equations 4 and 12, we can see the relationship between the kinetic parameters ($k^+$, $k^-$) and the bond graph parameters ($\kappa$, $K$):

By defining $\mathbf{Ln}$ as an element-wise logarithm operator, Equation 13 can be linearized and rewritten as the matrix equation:

In Pan (2019), the above equation was generalised using Equation 15.

where

$I_{n_r\times n_r}$ is an identity matrix of length $n_r$, while $n_r$ is the number of reactions. $N^{fT}$ and $N^{rT}$ are the transpose of forward and reverse stoichiometric matrices $N^{f}$ and $N^{r}$, respectively. The vectors of the forward and reverse kinetic rate constants $\mathbf{k^+}$, $\mathbf{k^-}$, and the vector of known constraints $K^c$ between the species defined in the matrix $N^c$ for SLC2A2 are shown in Equation 17.

The orders of the elements in $\mathbf{k^+}$ and $\mathbf{k^-}$ are the same order of the reactions, organized as columns in the matrices $N^f$ and $N^r$, shown in Table 2 and Table 3, while $N^c$ is organized by [number of species]$\times$[number of $K^c$]. Note that we do not need to add constraints in this case therefore both $K^c$ and $N^c$ are empty.

The diagonal matrix $\mathbf{W}$ (Equation 18) accounts for the volumes of compartments and the size is 10 ([number of reactions]$+$[number of species]). The typical blood cell volume $V_i=0.09$ ($pL$) according to McLaren et al. (1987) and we set the extracellular volume $V_o=0.09$ ($pL$) as well. Since the protein SLC2A2 does not exist in a compartment, we set volumes of corresponding conformations of the protein $V_1=V_2=V_3=V_4=1$ ($pL$) Pan, 2019. That is, their thermodynamic constants are not related to the volumes.

Given the above vectors and matrices, we can obtain bond graph parameters $\mathbf{\lambda}$ by matrix inversion (Equation 19),

where $\mathbf{M^+}$ is the Moore-Penrose pseudoinverse of $\mathbf{M}$ and $\mathbf{Exp}$ is the element-wise exponential. The resulting bond graph parameters based on Equation 19 are shown in Table 4.

3.2Parameters for steady state¶

Hunter et al. (2025) provided the analytical expression in Equation 20 Hunter et al., 2025 to calculate the parameters for the steady-state flux in Equation 19 Hunter et al., 2025 using the bond graph parameters in Table 4. We refer the readers to the primary paper for the analytical derivation process and calculation, while we explain here how we obtained the parameters from the steady-state data Lowe & Walmsley, 1986.

The Michaelis-Menten formulation of zero trans influx Lowe & Walmsley, 1986 of the transporter, i.e., set the intracellular concentration to be 0 ($mM$), is shown in Equation 20,

where $[A]_o$ is the extracellular concentration of glucose, the maximum flux $V_{oi}^{max}$ is calculated using Equation 21 with the concentration of glucose carrier molecules $[C]$ in human red blood cells of 6.67 ($\mu M$) Lowe & Walmsley, 1986.

The Michaelis-Menten constant of Equation 20 is calculated using Equation 22 Lowe & Walmsley, 1986, where $\frac{b}{a}$ is the dissociation constant of reaction $Re_3$.

When the intracellular molar amount of glucose is zero, the steady-state expression (Equation 19 in Hunter et al. (2025)) can be rearranged to Equation 23.

Substitute $q_o^A$ in the above equation with $[A]_oV_o$, and add the $V_E=1$ ($pL$) term to convert the unit from $fmol.s^{-1}$ to $mM.s^{-1}$ and rearrange it to Equation 24.

We obtain the following relationships in Equations 25 and 21 by comparing Equations 20 and 24.

Hence, we obtained the parameters $k_m^1$ and $v_m$ using Equations 27 and 28.

The Michaelis-Menten of zero trans efflux Lowe & Walmsley, 1986 i.e., set the extracellular concentration to be 0 ($mM$), is shown in Equation 29.

, where $[A_i]$ is the intracellular concentration of glucose, and the maximum flux $V_{io}^{max}$ is calculated using Equation 30 Lowe & Walmsley, 1986.

The Michaelis-Menten constant is calculated using Equation 31 Lowe & Walmsley, 1986, where $\frac{e}{f}$ is the dissociation constant of reaction $Re_4$.

When the extracellular molar amount of glucose is zero, the steady-state expression (Equation 19 in Hunter et al. (2025)) can be rearranged to Equation 32.

Substitute $q_i^A$ in the above equation with $[A]_iV_i$, and add the $V_E=1$ ($pL$) term to convert the unit from $fmol.s^{-1}$ to $mM.s^{-1}$ and rearrange it to Equation 33.

By comparing Equations 29 and 33, we obtained Equation 34.

Then we can calculate the parameter $k_m^2$ using Equation 35.

$k_m^3$ is calculated using Equation 20 in Hunter et al. (2025). The parameters are summarized in Table 5.

3.3Simulation results¶

We encoded the full bond graph model Hunter et al., 2025 in GLUT2_BG.cellml and the parameters in params_BG.cellml. To simulate the inward flux, we set the molar amount of intracellular glucose $q_{Ai}$ to be a very small value $0.09e-4$ ($fmol$), and varied the molar amount of extracellular glucose $q_{Ao}$ from $0.09e-4$ ($fmol$) to 2.25 ($fmol$). For each extracellular glucose value, we simulate the model for 250 seconds to get the steady-state flow rate. Similarly to simulate the outward flux, we set the molar amount of extracellular glucose $q_{Ao}$ to be a very small value $0.09e-4$ ($fmol$), and varied the molar amount of intracellular glucose $q_{Ai}$ from $0.09e-4$ ($fmol$) to 2.25 ($fmol$). For each intracellular glucose value, we simulate the model for 250 seconds to get the steady-state flow rate.

We encoded the Equations 20 and 29 in GLUT2_kinetic.cellml and simulated the zero trans influx and efflux by varying the extracellular glucose concentration and intracellular glucose concentration respectively from $1e-8$ ($mM$) to 25 ($mM$).

The steady-state model in Equations 19 and 20 in Hunter et al. (2025) were encoded in
GLUT2_ss_oi.cellml and GLUT2_ss_io.cellml. GLUT2_ss_oi.cellml sets the molar amount of intracellular glucose $q_{Ai}$ to be very small value $0.09e-8$ ($fmol$) and varies the molar amount of extracellular glucose $q_{Ao}$ from $0.09e-8$ ($fmol$) to 2.25 ($fmol$); GLUT2_ss_io.cellml sets the molar amount of extracellular glucose $q_{Ao}$ to be very small value $0.09e-8$ ($fmol$) and varies the molar amount of intracellular glucose $q_{Ai}$ from $0.09e-8$ ($fmol$) to 2.25 ($fmol$).

Figure 3 shows the steady-state fluxes from the full bond graph model and steady-state model in Hunter et al. (2025) compared to the zero trans influx (Equation 20) and efflux (Equation 29) in Lowe & Walmsley (1986). The plots in red used the parameters in Table 5 for Equation 19 in Hunter et al. (2025), while the magenta lines used the parameters in Table 4 to calculate the parameters for Equation 19 in Hunter et al. (2025) according to Equation 20 in Hunter et al. (2025).

We summarize the models, parameters, and corresponding simulation plots in Table 6. We have provided the Python scripts under the folder <src> to run the simulations and plot the data, while the SED-ML files in <Facilitated transporter\ CellMLV2> detail the simulation settings. To get the result in Figure 3, the Python scripts sim_GLUT2.py, mergeData_GLUT2.py, plot_GLUT2.py should run in sequence.

4SLC5A1 bond graph model parameterization¶

The SLC5A1 (SGLT1) uses the sodium gradient to drive glucose into the cell, typically when the transmembrane glucose gradient is insufficient to provide the required flux of glucose. The bond graph is shown in Figure 4. We parameterize the bond graph model to fit the data in Parent et al. (1992), and the kinetic model diagram is shown in Figure 5. Note that the notation and the parameter names in the kinetic diagram are different from the bond graph.

The kinetic parameters in in Parent et al. (1992) are listed in Table 7 and Table 8. The first column is the corresponding kinetic parameters for the reactions in the bond graph where the subscript is the reaction number. The original units in Parent et al. (1992) were $mole^{-2}.s^{-1}$, $mole^{-1}.s^{-1}$ or $s^{-1}$, while the units were changed to $M^{-2}.s^{-1}$, $M^{-1}.s^{-1}$ or $s^{-1}$ in Eskandari et al. (2005) where the model Parent et al., 1992 was reused. We found that using the units in Eskandari et al. (2005) gave the right dynamic outputs, so we used $M^{-2}.s^{-1}$, $M^{-1}.s^{-1}$ or $s^{-1}$ in this article and the primary paper Hunter et al., 2025.

4.1Bond graph parameters¶

We apply the same method Pan, 2019 to convert the kinetic parameters to bond graph parameters. Since we have described how to link kinetic parameters to bond graph parameters in Section 3.1, we directly solve Equation 15 to get the bond graph parameters for SLC5A1 without repeating the derivation process. First, we construct matrices $\mathbf{k^+}$, $\mathbf{k^-}$, $\mathbf{N^f}$, $\mathbf{N^r}$ and $\mathbf{W}$ to get Equation 15 for SLC5A1. The vectors of the forward and reverse kinetic rate constants $k^+$, $k^-$ are defined below and the values are from Table 7 and Table 8.

Note that we do not need to add constraints therefore both $K^c$ and $N^c$ are empty in this case. The forward and reverse stoichiometric matrices ($N^f,N^r$) are shown in Table 9 and Table 10, respectively. The first row lists the reactions while the first column denotes the species.

Based on McLaren et al. (1987) we set the volume of blood cells to $V_i=8.5 \times 10^{-2}$ ($pL$) and use the same value for the extracellular volume. The diagonal matrix $\mathbf{W}$ that accounts for the volumes of compartments is constructed in Equation 37.

Note that when preparing this Physiome manuscript we discovered a unit conversion error of the cell volume used to calculate parameter values in the Primary paper Hunter et al., 2025. This error scaled up the thermodynamical parameters of $Na$ and $Glc$ by 10⁷, i.e., $K_i^{Na}$ $K_o^{Na}$,$K_i^{Glc}$ and $K_o^{Glc}$ with corrected values given in Tables 11, 12 and 13. Since the product of the thermodynamical parameter and cell volume, e.g., $K_i^{Na}V_i$, determines the dynamics of the system, the scaling effect ($K_i^{Na}\cdot 10^7\cdot V_i\cdot 10^{-7}$) is canceled. Therefore, this change does not affect the model dynamics or simulation results.