#

# Ground state correlations on ground state densities and total binding
energies of ${}^{40}\mathrm{Ca}$, ${}^{48}\mathrm{Ca}$ and 208Pb

F. Minato<EMAIL_ADDRESS>Nuclear Data Center, Japan Atomic Energy
Agency, Tokai, Ibaraki 319-1195, Japan  RIKEN Nishina Center, Wako, Saitama
351-0198, Japan  H. Sagawa  RIKEN Nishina Center, Wako, Saitama 351-0198,
Japan  Center for Mathematical and Physics, the University of Aizu, Aizu-
Wakamatsu, Fukushima 965-8580, Japan  S. Yoshida  Science Research Center,
Hosei University, 2-17-1 Fujimi, Chiyoda, Tokyo 102-8160, Japan

###### Abstract

Neutron and proton densities of doubly-closed shell nuclei
${}^{40}\mathrm{Ca}$, ${}^{48}\mathrm{Ca}$ and 208Pb are studied based on a
Hartree-Fock model with SAMi, SAMi-J27 and SAMi-T energy density functionals
(EDFs). The ground state correlations (GSC) induced by isoscalar and isovector
phonons are also evaluated by the second order perturbation theory with a
self-consistent random phase approximation (RPA). We found that the interior
part of the ground state densities is reduced by the GSC in consistent with
the experimental data. On the other hand, the GSC enhances the neutron skin
thickness of ${}^{48}\mathrm{Ca}$ and 208Pb. The effect of GSC on the total
binding energy is also evaluated by the quasi-boson approximation. The effect
of the tensor interaction is found small on both the density distributions and
the binding energies.

ground state density, binding energy, Hartree-Fock approximation, ground state
correlation

##  1 Introduction

Mean field models have been quite successful to describe gross features of
ground state mean square radii and the total binding energies of many nuclei
in a wide region of the mass table Goriely2016 ; MinShi2019 ; Kortelainen2014
. Recent studies of the density distributions of doubly-closed shell nuclei
show appreciable differences between experimental and theoretical predictions
not only in the interior part but also the surface region Zenihiro ; Zenihiro2
. In the Hartree-Fock (HF) model for the doubly-closed shell nuclei, the
single-particle orbits are fully occupied up to the Fermi energy from the
bottom of the potential and no occupation probabilities of the orbits above
the Fermi energy. On the other hand, it has been reported by experiments
p-trans1 ; p-trans2  and theoretical large-scale shell model calculations in
40Ca Sagawa-PLB ; Brown ; Streets ; Shimizu ; Suzuki  that partial occupation
probabilities of the orbits take values between $0$ and $1$, especially around
the Fermi energy. These partial occupations are mostly induced by the many-
body correlations beyond the mean field models.

In this paper, we study the ground state correlations (GSCs) induced by
phonons calculated by a self-consistent Hartree-Fock (HF)+random phase
approximation (RPA) with Skyrme energy density functionals (EDFs). The GSC are
evaluated in the second-order perturbation model for the ground state
densities, and the quasi-boson approximation for the total binding energies of
${}^{40}\mathrm{Ca}$, ${}^{48}\mathrm{Ca}$ and 208Pb. A similar study has been
performed in Ref. Faessler1976 ; Waroquier1983  for 16O, 40,48Ca, and 208Pb.
This work, however, studies the GSC with a modern version of Skyrme
parameters, SAMi Roca-Maza2012 , SAMi-J27 Roca-Maza2013 , and SAMi-T Shen2019
EDFs and also focuses on the effect of the tensor force.

The paper is organized as follows. Section II is devoted to describe the
theoretical model for GSC. Numerical results are given in Section III. Summary
and future perspectives are given in Section IV.

##  2 Model for GSC

The ground state density distributions of doubly-closed shell nuclei are
commonly calculated by the Skyrme Hartree-Fock (SHF) model assuming spherical
symmetry. The density distribution of radial component is given as

$\rho_{q}(r)=\sum_{k=m,i}n_{k}\frac{2j_{k}+1}{4\pi}\left|\varphi_{k}(r)\right|^{2}\quad(q=p,\,n),$
(1)

where $n_{k}$, $\varphi_{k}$ and $j_{k}$ are the occupation probability, the
single-particle wave function of radial part and the total angular momentum
for state $k$, respectively. Here, we denote $m$ for particle states and $i$
for hole states. The occupation probabilities are $n_{m}=0$ and $n_{i}=1$ for
the HF density distribution. In case of the perturbed density distribution,
the occupation probabilities are modified within the number-operator method
Rowe1968c  to be

$\begin{split}n_{m}&=\frac{1}{2(2j_{m}+1)}\sum_{i\lambda
J\pi}(2J+1)\left|Y_{mi}(\lambda J^{\pi})\right|^{2}\\\
n_{i}&=1-\frac{1}{2(2j_{i}+1)}\sum_{m\lambda J\pi}(2J+1)\left|Y_{mi}(\lambda
J^{\pi})\right|^{2}.\end{split}$ (2)

where $Y_{mi}(\lambda J^{\pi})$ is the backward amplitude of the phonon
operator of random phase approximation. The label $\lambda$ and $J^{\pi}$
represent the $\lambda$-th excited RPA phonon state and its spin-parity,
respectively. The RPA phonon creation operator is defined as

$Q^{\dagger}(\lambda J^{\pi})=\sum_{mi}X_{mi}(\lambda
J^{\pi})a_{m}^{\dagger}a_{i}-Y_{mi}(\lambda J^{\pi})a_{i}^{\dagger}a_{m}.$ (3)

The coefficients of $X_{mi}(\lambda J^{\pi})$ and $Y_{mi}(\lambda J^{\pi})$
are determined by solving the RPA equation RingSchuck .

The rms radii are calculated by

$r_{q}\equiv\sqrt{\langle
r^{2}\rangle_{q}}=\sqrt{\frac{1}{N_{q}}4\pi\int\rho_{q}(r)r^{4}\,dr},$ (4)

where $N_{p}=Z$ and $N_{n}=N$. The neutron skin thickness is estimated by
$\Delta r_{np}=r_{n}-r_{p}$.

To calculate the unperturbed and perturbed density distributions numerically,
we adopt the SHF model in the coordinate space. To solve the RPA equation, we
take into account the single particle states up to $60$ MeV in energy, and
restrict unperturbed $1p$-$1h$ excitations energy up to $100$ MeV in the
model. The continuum states are discretized by a finite box with a size of
$R=20$ fm. We consider the spin-parity states up to $J^{\pi}=5^{\pm}$. We
confirm that the model space is enough to make the calculated occupation
probabilities of Eq. (2) converge. The residual interactions, the second-
derivative of SHF EDF with respect to densities, are self-consistently taken
into account in the RPA except for spin-orbit component, which is however not
significant. Within this model space, the energy-weighted sum rules are
satisfied over $99.5$% for all $J^{\pi}$.

##  3 Results

###  3.1 Density distributions and rms radii

Table 1: Proton and neutron rms radii of ${}^{40}\mathrm{Ca}$, ${}^{48}\mathrm{Ca}$, and 208Pb. The neutron skin is given as $\Delta r_{np}\equiv r_{n}-r_{p}$. HF denotes the results of HF calculations, while HF+GSC includes also the ground state correlations. The experimental data are taken from Zenihiro  for ${}^{40}\mathrm{Ca}$, ${}^{48}\mathrm{Ca}$ and from Zenihiro2  for 208Pb. The experimental errors of $\Delta r_{np}$ include both statistical errors and model uncertainties. The radii are given in unit of fm. |  | SAMi | SAMi-J27 | SAMi-T |
---|---|---|---|---|---
Nuclei |  | HF | HF+GSC | HF | HF+GSC | HF | HF+GSC | exp.
| $r_{n}$ | 3.343 | 3.773 | 3.344 | 3.810 | 3.339 | 3.808 | 3.375
40Ca | $r_{p}$ | 3.390 | 3.853 | 3.391 | 3.870 | 3.386 | 3.845 | 3.385
| $\Delta r_{np}$ |  $-$0.047 |  $-$0.080 |  $-$0.047 |  $-$0.060 |  $-$0.046 |  $-$0.037 | $-0.010^{+0.049}_{-0.048}$
| $r_{n}$ | 3.612 | 4.003 | 3.588 | 4.006 | 3.589 | 4.004 | 3.555
48Ca | $r_{p}$ | 3.436 | 3.770 | 3.444 | 3.801 | 3.424 | 3.781 | 3.387
| $\Delta r_{np}$ | 0.176 | 0.233 | 0.144 | 0.205 | 0.166 | 0.223 | $0.168^{+0.052}_{-0.055}$
| $r_{n}$ | 5.610 | 5.750 | 5.580 | 5.731 | 5.574 | 5.730 | 5.653
208Pb | $r_{p}$ | 5.463 | 5.555 | 5.456 | 5.553 | 5.421 | 5.517 | 5.442
| $\Delta r_{np}$ | 0.147 | 0.195 | 0.123 | 0.178 | 0.153 | 0.213 | $0.211^{+0.054}_{-0.063}$

Table 1 shows nuclear proton and neutron rms radii of ${}^{40}\mathrm{Ca}$,
${}^{48}\mathrm{Ca}$ and 208Pb calculated by using a modern version of Skyrme
parameters, SAMi, SAMi$-$J27, and SAMi-T EDFs. Both the results of HF and
HF+GSC calculations are listed and compared with experimental data. The tensor
interactions are included in SAMi-T EDF. The proton radii or equivalently the
charge radii of ${}^{40}\mathrm{Ca}$ and 208Pb are in general included in the
optimized data set. Therefore it is not a surprise that HF results reproduce
well the proton radii of these nuclei. The neutron radii of
${}^{40}\mathrm{Ca}$ and 208Pb are slightly smaller than the experimental
ones. Including the GSC, the neutron skin of 208Pb increases and becomes
closer to the experimental value compared with the HF result. On the other
hand, the neutron skin of ${}^{40}\mathrm{Ca}$ becomes a larger negative value
with the GSC in the cases of SAMi and SAMi-J27, while SAMi-T with the tensor
term gives slightly better agreement with the experimental value $\Delta
r_{np}=-0.01$ fm. In ${}^{48}\mathrm{Ca}$, all the three HF results in Table 1
give slightly larger radii than the experimental one for both protons and
neutrons, while the neutron skin is well reproduced by the HF calculations,
especially in the case of SAMi-T. The GSC enlarges both the proton and the
neutron radii and the quantitative agreement becomes worse than HF results by
about 0.4 fm. Experimental data of the neutron skin of 208Pb is still
controversial. The dipole polarizability experiment gives $\Delta
r_{np}=0.156^{+0.025}_{-0.021}$ fm Tamii , while the parity violation electron
scattering experiment PREXII were analyzed by two different models and the
neutron skin was extracted to be $\Delta r_{np}=0.283\pm 0.071$ fm PREXII  and
$\Delta r_{np}=0.19\pm 0.02$ fm PG-PREX , respectively. The present
calculations, both HF and HF+GSC, are unlikely to reproduce the experimental
data of PREXII reported in ref. PREXII .

Figure 1: (color online) Neutron (left panels) and proton (right panels)
density distributions of ${}^{40}\mathrm{Ca}$. Upper panels show the HF
results with three SAMi EDFs, SAMi, SAMi-J27 and SAMi-T, respectively. The
SAMi-T has the tensor terms, while the SAMi and SAMi-J27 have no tensor terms.
The black solid lines show experimental data taken from Ref. electron-sca  for
protons and from Ref. Zenihiro  for neutrons. The shaded area of experimental
neutron density shows experimental uncertainties of both statistical and
systematic errors. The neutron density distribution denoted SAMi-J27mod in the
upper left panel is obtained by optimizing the occupation probabilities of
neutrons as given in Table 2.  Figure 2: (color online) Neutron and proton
density distributions of ${}^{40}\mathrm{Ca}$ in the log scale. The upper
panels show the HF results, while the lower panels are the HF+GCS results.
Experimental data are shown by a black solid curve. The model dependence of
EDFs is not visible in this scale.

Figure 1 shows comparisons between calculated and experimental neutron (left
panel) and proton (right panel) densities of 40Ca. Upper panels show the HF
results, and the lower panels correspond to the results of HF+GSC. The HF
densities of both neutrons and protons somewhat overestimate in the interior
part. The EDF model dependence is rather small: the SAMi and SAMi-J27 give
almost identical results for the density distributions, while SAMi-T enhances
the interior region of density more than other two interactions. On the
contrarily, the calculated results slightly underestimate the shoulder part at
around $3$ fm. With the GSC, the central part is hindered and become closer to
the experimental density distribution, especially in the case of proton
density. On top of that, the shoulder part is also quenched and small amount
of density is shifted to $r>6$ fm region. In the upper left panel of Fig. 1,
the neutron density distribution denoted SAMi-J27mod is shown by optimizing
the occupation probabilities of neutrons as given in Table 2 in order to
reproduce the experimental density. The difference between HF and HF+GSC at
$r>6.5$ fm is clearly seen in Fig. 2 where the densities are plotted in log
scale. In this scale, the interaction dependence does not appear in the entire
region of coordinate space.

Figure 3 shows the neutron density distribution of 40Ca without and with GSC
multiplied by a factor $4\pi r^{2}$. The differences between HF and HF+GSC
appear in the central region $r<2$ fm, the peak height at $r=3.2$ fm, in the
tail region 5 fm$<r<$6 fm, and also at the very low density region $r>6.5$ fm.
We can see small improvements of agreements by GSC at the central region and
also at the surface. On the other hand, the peak height is quenched by GSC and
a small amount of density is pushed to outside of the surface region in
contrast to the experimental data.

Figure 3: (color online) Neutron density distributions of ${}^{40}\mathrm{Ca}$
with a factor of $4\pi r^{2}$. The left panel shows the HF results, while the
right panel are the HF+GCS results. Experimental data are shown by a black
solid curve. The model dependence of EDFs is not visible in this scale. Figure
4: (color online) The same as Fig. 1, but for ${}^{48}\mathrm{Ca}$.

Figure 4 shows comparisons between calculated and experimental densities of
neutrons (left panel) and protons (right panel) of 48Ca. The essential
features are the same as those of 40Ca. One peculiar difference is the
shoulder peak of neutron density due to the $f_{7/2}$ neutron occupation in
48Ca. The HF densities of both neutrons and protons slightly overestimate in
the interior part. The EDF model dependence is again rather small for 48Ca.
The calculated results underestimate the shoulder part at around $3$ fm. With
the GSC, the central part is hindered and become closer to the experimental
density distribution, or even underestimate in the case of proton density. The
shoulder part is also quenched and small amount of density is shifted to $r>6$
fm region.

Figure 5: (color online) The same as Fig. 1, but for ${}^{208}\mathrm{Pb}$.
Experimental data taken from Ref. electron-sca  for protons and from Ref.
Zenihiro2  for neutrons.

The densities of 208Pb are shown in Fig. 5. For the proton density (right
panel), the HF results reproduce well the plateau between $r=2.0$ fm and
$r=5.5$ fm, while they show larger values at the central part at $r<1.5$ fm
than the experimental data. The GSC effect on the proton density of 208Pb is
relatively small even compared with 40Ca and 48Ca. Nevertheless, we can see
that the HF+GSC improves slightly the agreement at the interior part of proton
density at the right panels of Fig. 5. For the neutron density, the HF results
give slightly larger density for the plateau between $r=1.0$ and $r=5.5$ fm,
while underestimate the central part $r<1$ fm. These disagreements are
somewhat improved in the case of HF+GSC results and the calculated results are
within the experimental uncertainties at $r<5$ fm. he density distribution at
$r>8$ fm becomes larger than that of HF similar to the case of 40Ca shown in
Fig. 2.

###  3.2 Particle occupation numbers of $sd$-shell and $pf$-shell orbits in
40Ca

In Table 2, the occupation probabilities of various calculations are listed
with the empirical proton and neutron occupation probabilities in 40Ca. Skyrme
HF-Bogolyubov calculations are performed with SAMi interaction and the volume-
type and mixed-type pairing interactions with the code HFBTHO Sagawa-PLB . We
found that there are essentially no occupation probabilities $v_{j}^{2}$ in
$pf$-shell orbits because of a large shell gap between $1d_{3/2}$ and
$1f_{7/2}$ orbits. The GSCs with SAMi and SAMi-T give small amounts of the
occupation probabilities in the $pf$-shell orbits,
$\sum_{pf}(2j+1)v_{j}^{2}\approx$ 0.3. There is no essential difference
between SAMi and SAMi-T without and with the tensor terms.

The results of several Large Scale Shell Model (LSSM) calculations are also
listed in Table 2, which depend largely on the model space and the adopted
interactions. In Refs. Brown ; Shimizu ; Suzuki , LSSM calculations have been
performed including many-particle many-hole excitations from $sd$-shell to
$pf$-shell configurations. In Ref. Brown , the active model space is
$(1d_{3/2},1f_{7/2},2p_{3/2})$, while the $1d_{3/2},2s_{1/2}$ orbits are also
included in Ref. Shimizu . Shell model calculations with the configurations of
2-particle 2-hole (2p-2h) and 4p-4h excitations from the closed shell core of
40Ca are also performed in Ref. Suzuki . The full $sd$ and $pf$ shell orbits
are involved with $sdpf$-mu shell model effective interactions. The occupation
probabilities of $1f$-orbits go up to $\sum_{f}(2j+1)v_{j}^{2}\approx 0.6$,
but those of $p$-shell orbits are rather small as less than 0.1. The empirical
occupation numbers are also listed obtained from proton transfer reactions
p-trans1 , and also from the analysis of neutron density distributions denoted
as SAMi-J27mod in Fig. 1 Sagawa-PLB . The proton transfer data suggest the
same amount of the proton occupation probabilities in $pf$-shell orbits as the
shell model predictions. On the other hand, the experimental neutron density
distribution requires the occupation probability of $f$-orbits to be
$\sim$1.2, which is two times larger than the proton transfer experiments.

Table 2: Particle occupation numbers, $(2j+1)v_{j}^{2}$, for neutron $sd$-shell and $pf$-shell configurations of 40Ca. The column with a bar (—) is not involved in the shell model calculations or analysis of experimental data. Experimental data of proton transfer reaction are taken from Ref. p-trans1 . Model | $1s_{1/2}$ | $2p_{3/2}$ | $2p_{1/2}$ | $1d_{5/2}$ | $2s_{1/2}$ | $1d_{3/2}$ | $1f_{7/2}$ | $2p_{3/2}$ | $2p_{1/2}$ | $1f_{5/2}$
---|---|---|---|---|---|---|---|---|---|---
HFB (SAMi) | 2.00 | 4.00 | 2.00 | 6.00 | 2.00 | 4.00 | 0.00 | 0.00 | 0.00 | 0.00
HF+GSC(SAMi) | 1.964 | 3.908 | 1.949 | 5.761 | 1.843 | 3.792 | 0.206 | 0.049 | 0.010 | 0.040
HF+GSC(SAMi-T) | 1.968 | 3.905 | 1.954 | 5.753 | 1.845 | 3.776 | 0.201 | 0.052 | 0.009 | 0.076
GSC (RPA) Lenske  | 1.92 | 3.76 | 1.88 | 5.58 | 1.8 | 3.7 | 0.24 | 0.16 | 0.08 | —
$dpf$-shell Brown  | — | — | — | — | — | 3.30 | 0.63 | 0.07 | — | —
$sdpf$-msd4 Shimizu  | — | — | — | 5.902 | 1.908 | 3.477 | 0.617 | 0.096 | — | —
$sdpf$-mu $(2p$-$2h)$ Suzuki  | — | — | — | 5.864 | 1.933 | 3.845 | 0.191 | 0.040 | 0.020 | 0.107
$sdpf$-mu $(4p$-$4h)$ Suzuki  | — | — | — | 5.727 | 1.854 | 3.660 | 0.421 | 0.089 | 0.041 | 0.208
exp. (p-transfer) p-trans1  | — | — | — | 6.0 | 1.70 | 3.59 | 0.56 | 0.15 |  |
SAMi-J27mod Sagawa-PLB  | 2.00 | 4.00 | 2.00 | 6.00 | 1.342 | 3.431 | 1.227 |  |  |

###  3.3 Correlation energies for the total binding energies

The correlation energies due to the ground state correlations are calculated
within the framework of RPA and quasi-boson approximation Rowe1968c :

$E_{\rm corr}=\sum_{\lambda J\pi}(2J+1)E(\lambda
J^{\pi})\sum_{mi}\left|Y_{mi}(\lambda J^{\pi})\right|^{2}.$ (5)

The sum of phonon states $\lambda$ in Eq. (5) is restricted to those absorb a
fraction of the non-energy-weighted sum rule larger than $1$% for isoscalar or
isovector transition strength. The results are tabulated in Table 3. As is
expected, the HF calculations give already reasonable agreement with the
experimental total binding energies since all the parameters are optimized to
reproduce a set of experimental data in which the total binding energies of
the doubly-closed shell nuclei are included. The correlation energies give
about $1$% additional energies to the binding energies of 40Ca and 48Ca, while
it is much smaller in 208Pb to be about $0.2$%. As seen in the table, the
correlation energies do not give further improvement. For this reason,
together with the rms radii, one needs to readjust the Skyrme parameters
including the correlation effects for the consistent description of the ground
state properties in the beyond mean field model.

Table 3: The HF total binding energies and correlation energies due the ground state correlations of ${}^{40,48}\mathrm{Ca}$ and 208Pb. The experimental data are taken from AME2020 AME2020 .The correlation energies are calculated by using Eq. (5). The value is given by a unit of MeV. | SAMi | SAMi-J27 | SAMi-T |
---|---|---|---|---
Nuclei | Etot | Ecorr | Etot | Ecorr | Etot | Ecorr | exp.
40Ca | 347.02 | 4.25 | 344.39 | 4.10 | 343.68 | 4.78 | 342.05
48Ca | 415.46 | 4.93 | 415.45 | 5.89 | 416.73 | 4.76 | 416.00
208Pb | 1635.86 | 3.95 | 1635.54 | 3.05 | 1635.39 | 3.41 | 1636.43

##  4 Summary and future perspectives

We studied the effects of GSC on the densities and the total binding energies
of doubly-closed shell nuclei 40Ca, 48Ca, and 208Pb with modern Skyrme EDFs,
SAMi, SAMi-J27 and SAMi-T. For the ground state densities, the GSC reduces
slightly in the interior part of HF densities of both neutrons and protons,
and gives consistent results in comparisons with the experimental data of 40Ca
and 48Ca. The GSC decreases also the shoulder part of density at around $r=3$
fm, and removes a small amount of density to the region $r>6$ fm. The GSC
effect is smaller in 208Pb and quantitatively the same as 40Ca, and 48Ca.

The GSC on the total binding energy is also discussed within the quasi-boson
approximation, and found to increase about $1$% for 40,48Ca and $0.2$% for
208Pb of the total binding energy. In order to describe these ground state
observables with beyond mean field models, we need to refit Skyrme EDFs
including the GSC effects for the consistent descriptions.

## Acknowledgement

FM and HS thank YIPQS long-term workshop “Mean-field and Cluster Dynamics in
Nuclear Systems 2022 (MCD2022)" May 9 – June 17, 2022 held at Yukawa
Institute, Kyoto University, where the collaboration was started. We thank J.
Zenihiro and T. Uesaka for providing the experimental data. We thank also T.
Naito for useful discussions on the HF calculations. This work is supported by
JSPS KAKENHI Grant Number JP19K03858.JSPS.

## References

  * (1) S. Goriely, N. Chamel, and J.M. Pearson, Phys. Rev. C 93, 034337 (2016).
  * (2) Min Shi, Zhong-Ming Niu, and Hao-Zhao Liang, Chin. Phys. C 43, 074104 (2019).
  * (3) M. Kortelainen, J. McDonnell, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sarich, N. Schunck, S.M. Wild, D. Davesne, J. Erler, and A. Pastore, Phys. Rev. C 89, 054314 (2014).
  * (4) J. Zenihiro, T. Uesaka, S. Yoshida and H. Sagawa, Prog, Theo. Exp. Phys. 2021, 023D05 (2021)
J. Zenihiro et al., arXiv:1810.11796 (2018).

  * (5) J. Zenihiro et al., Phys. Rev. C82, 044611 (2010).
  * (6) P. Doll, G.J. Wagner, K.T. Knöpfle and G. Mairle, Nucl. Phys. A 263, 210 (1976).
  * (7) F. Malaguti, A. Uguzzoni, E. Verondini and P. E. Hodgson, Nucl. Phys. A 297, 287 (1978); ibid., Nuovo Cim. A 49, 412 (1979).
  * (8) B. A. Brown, S. E. Massen, and P. E. Hodgson, J. Phys. G: Nucl. Phys. 5, 1655 (1979).
  * (9) J. Streets, B.A. Brown and P.E. Hodgson, J. Phys. G: Nucl. Phys. 8, 839 (1982).
  * (10) N. Shimizu, Y. Utsuno, T. Ichikawa et al., private communications (2022).
  * (11) Toshio Suzuki, private communications (2022).
  * (12) H. Sagawa, S. Yoshida, T. Naito, T. Uesaka, J. Zenihiro, J. Tanaka, Physics Letters B829 137072 (2022).
  * (13) A. Faessler, S. Krewald, A. Plastino, and J. Speth, Z. Physik A 276, 91 (1976).
  * (14) M. Waroquier, J. Bloch, G. Wenes, and K. Heyde, Phys. Rev. C 28, 1791 (1983).
  * (15) X. Roca-Maza, G. Colò, and H. Sagawa, Phys. Rev. C 86, 031306(R) (2012).
  * (16) X. Roca-Maza, M. Brenna, B. K. Agrawal, P. F. Bortignon, G. Colò, Li-Gang Cao, N. Paar, and D. Vretenar, Phys. Rev. C 87, 034301 (2013).
  * (17) S. Shen, G. Coló, and X. Roca-Maza, Phys. Rev. C 99, 034322 (2019).
  * (18) D.J. Rowe, Phys. Rev. 175, 1283 (1968).
  * (19) P. Ring and P. Schuck, The Nuclear Many-Body Problem, Springer-Verlag, New York, 1980.
  * (20) H. Lenske and J. Wambach, Phys. Lett. B 249, 377 (1990).
  * (21) A. Tamii et al., Phys. Rev. Lett. 107, 062502 (2011).
  * (22) D. Adhikari et al. (PREX Collaboration), Phys. Rev. Lett. 126, 172502 (2021).
  * (23) P.-G. Reinhard, X. Roca-Maza, and W. Nazarewicz Phys. Rev. Lett. 127, 232501 (2021).
  * (24) H. de Vries, C. W. de Jager, and C. de Vries, At. Data Nucl. Data Tables 36, 495 (1987).
  * (25) W. J. Huang, M. Wang, F.G. Kondev, G. Audi, and S. Naimi, Chin. Phys. C45, 030002 (2021); Chin. Phys. C45, 030003 (2021).


# Linear Response Theory of Evolved Metabolic Systems

Jumpei F. Yamagishi  Tetsuhiro S. Hatakeyama  Department of Basic Science, The
University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan

###### Abstract

Predicting cellular metabolic states is a central problem in biophysics.
Conventional approaches, however, sensitively depend on the microscopic
details of individual metabolic systems. In this Letter, we derived a
universal linear relationship between the metabolic responses against nutrient
conditions and metabolic inhibition, with the aid of a microeconomic theory.
The relationship holds in arbitrary metabolic systems as long as the law of
mass conservation stands, as supported by extensive numerical calculations. It
offers quantitative predictions without prior knowledge of systems.

Metabolism is the physicochemical basis of life. Understanding its behavior
has been a major goal of biophysics ilker2019modeling ; kar2003collapse ;
yamagishi2020advantage ; zwicker2017growth . At the same time, the prediction
of cellular metabolic states is a central problem in biology. In particular,
prediction of the responses of metabolic systems against environmental
variations or experimental operations is essential for manipulating metabolic
systems to the desired states in both life sciences and in applications such
as matter production in metabolic engineering stephanopoulos1998metabolic  and
the development of drugs targeting cellular metabolism murima2014targeting ;
drug-induced-hyperthermia ; martinez2017cancer .

Previous studies have mainly attempted to predict the metabolic responses by
predicting the metabolic states before and after perturbations, and they
require building an ad hoc model for each specific metabolic system. In
systems biology, constraint-based modeling (CBM) has often been used to
predict the cellular metabolic states palsson2015systems ; klipp2016systems ;
warren2007duality . In this method, the intracellular metabolic state is
predicted by solving an optimization problem of models of metabolic systems,
including a detailed description of each metabolic reaction. To construct the
optimization problem, metabolic systems of cells are assumed to be optimized
through (sometimes artificial) evolution for some objectives
heinrich1998modelling ; klipp2016systems ; palsson2015systems , e.g.,
maximization of the growth rate in reproducing cells such as cancer cells and
microbes OMbook  and maximization of the production of some molecules in
metabolically engineered cells portnoy2011adaptive . Indeed, metabolic systems
of reproducing cells exhibit certain ubiquitous phenomena across various
species, and those phenomena can be explained as a result of optimization
under physicochemical constraints OMbook . Although the assumption of optimal
metabolic regulation seems acceptable, knowing the true objective function of
cells, which is essential to making a model for CBM, remains nearly
impossible. Besides, even with remarkable progress in omics research, fully
reconstructing metabolic network models for each individual species or cell of
interest is still a challenge. Moreover, the numerical predictions are
sensitive to the details of the concerned constraints and the objective
functions selected bonarius1997flux ; raman2009flux ; schnitzer2022choice ;
Heinemann-Gibbs2019 . Therefore, new methods independent of the details of
metabolic systems are required.

Instead of metabolic states themselves, here, we focus on the responses of
metabolic systems to perturbations. At first glance, such prediction is
seemingly more difficult than predicting the cellular metabolic states because
it seems to require information not only on the steady states but also on
their neighborhoods. However, from another perspective, to predict only the
metabolic responses, we may need to understand the structure of only a limited
part of the state space of feasible metabolic states. In contrast, we must
seek the whole space to predict the metabolic states themselves. If
optimization through evolution and some physicochemical features unique to
metabolic systems constrain the behavior in the state space, there might be
universal features in the responses of metabolic systems to perturbations,
independent of system details, as in the linear response theory in statistical
mechanics onsager1931reciprocal ; kubo1957statistical ; green1954markoff .

In this Letter, we demonstrate a universal property of intracellular metabolic
responses in the optimized metabolic regulation, using a microeconomic theory
Varian ; Lancaster1966 ; yamagishi2021microeconomics . By introducing a
microeconomics-inspired formulation of metabolic systems, we can take
advantage of tools and ideas from microeconomics such as the Slutsky equation
that describes how consumer demands change in response to income and price. We
thereby derive quantitative relations between the metabolic responses against
nutrient abundance and those against metabolic inhibitions, such as the
addition of metabolic inhibitors and leakage of intermediate metabolites; the
former is easy to measure in experiments while the latter may not be. The
relations universally hold independent of the details of metabolic systems as
long as the law of mass conservation holds. Our theory is applicable to any
metabolic system and will provide quantitative predictions on the
intracellular metabolic responses without detailed prior knowledge of
microscopic molecular mechanisms and cellular objective functions.

Microeconomic formulation of metabolic regulation.— We first provide a
microeconomic formulation of optimized metabolic regulation, which is
equivalent to linear programming problems in CBM (Fig. 1).

We denote the set of all chemical species (metabolites) and that of all
constraints by ${\mathcal{M}}$ and $\mathcal{C}$, respectively. $\mathcal{C}$
can reflect every type of constraints such as the allocation of proteins
scott2011bacterial ; OM , intracellular space vazquez2010catabolic , membrane
surfaces memRealEstate2 , and Gibbs energy dissipation Heinemann-Gibbs2019  as
well as the bounds of reaction fluxes.

In the microeconomic formulation, variables to be optimized are the fluxes of
metabolic pathways, whereas they are fluxes of reactions in usual CBM
approaches; a metabolic pathway is a linked series of reactions and thus
comprises multiple reactions. The sets of reactions and pathways are denoted
by ${\mathcal{R}}$ and ${\mathcal{P}}$, respectively. Let us then consider two
stoichiometry matrices for reactions and pathways, $S$ and $K$, respectively
(see also SM, Table 1). For chemical species $\alpha\,(\in{\mathcal{M}})$,
$|S_{\alpha i}|$ represents the number of units of species $\alpha$ produced
if $S_{\alpha i}\,{>}\,0$ and consumed if $S_{\alpha i}\,{<}\,0$ in reaction
$i$; whereas if $\alpha$ denotes a constraint ($\alpha\in{\mathcal{C}}$),
$S_{\alpha i}$ is usually negative and $|S_{\alpha i}|$ represents the number
of units of constraint $\alpha$ required for reaction $i$. The stoichiometry
matrix $K$ for metabolic pathways ${\mathcal{P}}$ is also defined similarly.
Throughout the Letter, we use indices with primes such as ${i^{\prime}}$ to
denote pathways and those without primes such as $i$ to denote reactions, and
$|S_{\alpha i}|$ and $|K_{\alpha{i^{\prime}}}|$ are called input (output)
stoichiometric coefficients of reaction $i$ and pathway ${i^{\prime}}$,
respectively, if $S_{\alpha i}$ and $K_{\alpha{i^{\prime}}}$ are negative
(positive).

Cells are assumed to maximize the flux of some objective reaction
$o\,(\in{\mathcal{R}})$ such as biomass synthesis in reproducing cells and
ethanol or ATP synthesis in metabolically engineered cells.

We define the set of the species consumed in and the components required for
reaction $o$ as objective components
${\mathcal{O}}\,(\subset{\mathcal{M}}\cup{\mathcal{C}})$, and thus $S_{\alpha
o}$ for each objective component $\alpha\,(\in{\mathcal{O}})$ is negative.
Because the reactants of a reaction cannot be compensated for each other due
to the law of mass conservation yamagishi2021microeconomics ; liao2020modeling
; roy2021unifying , the flux of objective reaction $o$, i.e., the objective
function, is limited by the minimum available amount of objective components
${\mathcal{O}}$ as follows:

$\displaystyle\Lambda({{\bf
f}}):=\min_{\alpha\in{\mathcal{O}}}\left[\frac{1}{-S_{\alpha
o}}\left(\sum_{j^{\prime}\in{\mathcal{P}}}K_{\alpha
j^{\prime}}f_{j^{\prime}}+I_{\alpha}\right)\right],$ (1)

where ${{\bf f}}\,{=}\,\\{f_{i^{\prime}}\\}_{{i^{\prime}}\in{\mathcal{P}}}$
represents the fluxes of metabolic pathways. The arguments of the above min
function represent biologically different quantities: if $\alpha$ is a species
($\alpha\in{\mathcal{M}}$), $I_{\alpha}$ is its intake flux and
$\sum_{j^{\prime}}K_{\alpha j^{\prime}}f_{j^{\prime}}$ represents its total
production rate, while if $\alpha$ is a constraint ($\alpha\in{\mathcal{C}}$),
$I_{\alpha}$ is the total capacity for constraint $\alpha$ and
$\sum_{j^{\prime}}K_{\alpha j^{\prime}}f_{j^{\prime}}+I_{\alpha}$ is the
amount of $\alpha$ that can be allocated to the objective reaction.

The optimized solution $\hat{{{\bf f}}}$ is determined as a function of $K$
and I with the following constraints for the available pathway fluxes ${{\bf
f}}$:

$\displaystyle-\sum_{j^{\prime}\in{\mathcal{P}}}K_{\alpha
j^{\prime}}f_{j^{\prime}}\leq
I_{\alpha}.\quad{(\alpha\in\mathcal{E}\cup{\mathcal{C}})}$ (2)

Here, $\mathcal{E}\,(\subset{\mathcal{M}})$ denotes the set of exchangeable
species that are transported through the cellular membrane. That is, the above
constraints reflect that the total consumption of species cannot exceed their
intakes. If species $\alpha$ is produced by objective reaction $o$, the intake
effectively increases and $S_{\alpha o}\Lambda$ is added to the right-hand
side of Eq. (2), although this is not the case for most species.

This optimization problem (1-2) can be interpreted as a microeconomic problem
in the theory of consumer choice Varian ; Lancaster1966 ;
yamagishi2021microeconomics , considering $\Lambda({{\bf f}})$ as the utility
function. By focusing on an arbitrary component $\nu$, one of inequalities (2)
serves as the budget constraint for $\nu$ if $K_{\nu j^{\prime}}\,{\leq}\,0$
for all pathways $j^{\prime}$, while the remaining inequalities in Eq. (2)
then determine the solution space [Fig. 1(a)]: for example, if we choose
glucose as $\nu$, the corresponding inequality in Eq. (2) represents carbon
allocation. Here, the maximal intake $I_{\nu}$ of $\nu$ corresponds to the
income, and the input stoichiometric coefficient for each pathway,
$p_{j^{\prime}}^{\nu}:\,{=}\,-K_{\nu j^{\prime}}$, serves as the price of
pathway $j^{\prime}$ in terms of $\nu$.

Figure 1:  Schematic illustration. (a) (left) Metabolic CBM formulation with
reaction fluxes ${{\bf v}}$ as variables. The solution subspace (convex set of
possible allocations), called the flux cone, is shown in pink. (right)
Microeconomic formulation with pathway fluxes ${{\bf f}}$ as variables and an
objective flux $\Lambda$. The pink area in ${{\bf f}}$-plane (bottom surface)
represents the solution subspace, whereas the blue plane vertical to ${{\bf
f}}$-plane is the budget constraint for a component $\nu$. The blue points
$\hat{{{\bf v}}}$ and $\hat{{{\bf f}}}$ represent the optimized fluxes of
reactions and pathways, respectively. Given ${{\bf v}}\,{=}\,P{{\bf f}}$ with
pathway matrix $P$, both formulations are equivalent optimization problems
(see Supplemental Material (SM), Sec. S1 for details and Sec. S2 and Fig. S1
for a simple example). (b) Liner relation between the metabolic responses
against changes in nutrient conditions (yellow) and those against metabolic
inhibitions (green) [Eq. (3)).

Relation between responses of pathway fluxes to nutrient abundance and
metabolic inhibition.— Because Eqs. (1-2) can be interpreted as a
microeconomic optimization problem, we can apply and generalize the Slutsky
equation in the theory of consumer choice Varian . The equation shows the
relationship between changes in the optimized demands for goods in response to
income and price. In metabolism, it corresponds to the relationship between
the responses of optimal pathway fluxes $\hat{{{\bf f}}}$ (see SM, Sec. S4 for
derivation):

$\displaystyle\frac{\partial\hat{f}_{i^{\prime}}(K,{\bf I})}{\partial
p^{\nu}_{j^{\prime}}}=-\hat{f}_{j^{\prime}}(K,{\bf
I})\frac{\partial\hat{f}_{i^{\prime}}(K,{\bf I})}{\partial I_{\nu}}.$ (3)

The right-hand side represents the responses of pathway ${i^{\prime}}$ against
increases in $I_{\nu}$, whereas the left-hand side represents those against
metabolic inhibitions in pathway ${j^{\prime}}$ because the metabolic price
$p^{\nu}_{j^{\prime}}\,{=}\,-K_{\nu{j^{\prime}}}$ quantifies the inefficiency
of conversion from substrate $\nu$ to endproducts in pathway ${j^{\prime}}$
yamagishi2021microeconomics .

Figure 2: Responses of the optimized pathway fluxes ${\bf\hat{f}}$. (a)
Responses to metabolic inhibitions, $\Delta\hat{f}_{i^{\prime}}(\Delta
K_{\mathrm{Glc},{j^{\prime}}})/\Delta p^{\mathrm{Glc}}_{j^{\prime}}$, are
plotted against the nutrient responses,
$\hat{f}_{j^{\prime}}\Delta\hat{f}_{i^{\prime}}(\Delta
I_{\mathrm{Glc}})/\Delta I_{\mathrm{Glc}}$. All different shapes and colors of
markers represent different ${i^{\prime}}$ and ${j^{\prime}}$, respectively.
$I_{\mathrm{Glc}}\,{=}\,5\,\mathrm{[mmol/gDW/h]}$. (b) $13$ active extreme
pathways, computed using efmtool efmtool , are shown. Colors correspond to
those of the markers for manipulated pathways ${j^{\prime}}$ in panel (a). The
whole metabolic network of the E. coli core model is shown in gray.

The derivation of Eq. (3) relies solely on the law of mass conservation, i.e.,
the reactants of a reaction cannot be compensated for each other. Because the
law of mass conservation stands in every chemical reaction, the relation (3)
of the two measurable quantities must hold in arbitrary metabolic systems as
long as their metabolic regulation is optimized for a certain objective. In
particular, the case ${i^{\prime}}\,{=}\,{j^{\prime}}$ will be useful: it
indicates that measuring the responses of a pathway flux to changes in the
nutrient environment provides quantitative predictions of the pathway’s
responses to metabolic inhibition or activation, and vice versa.

To confirm the validity of Eq. (3), we numerically solved the optimization
problems (1-2) with pathway fluxes ${{\bf f}}$ as variables using the E.coli
core model orth2010reconstruction ; palsson2015systems  and randomly chosen
stoichiometric coefficients for the single constraint (Fig. 2). In this
numerical calculation, metabolic pathways from exchangeable species to
objective components are chosen as linear combinations of extreme pathways or
elementary flux modes schilling2000theory  for stoichiometry without objective
reaction $o$ [Fig. 2(b)], although the above arguments do not depend on the
specific choices of metabolic pathways (see SM, Sec. S3 for details). As shown
in Fig. 2(a), the linear relation (3) between metabolic responses is indeed
satisfied. Notably, it is satisfied regardless of the number and type of
constraint(s) ${\mathcal{C}}$, whereas the metabolic states themselves can
sensitively depend on the concerned constraints and environmental conditions.

Relation between responses of reaction fluxes.— Although Eq. (3) generally
holds for arbitrary metabolic pathways, it may be experimentally easier to
manipulate a single metabolic reaction. Manipulation of a single reaction can
affect multiple pathways because they are often tangled via a common reaction
in the metabolic network. Thus, we should consider the contributions of
multiple pathways. The simplest way for this is to sum up Eq. (3) for all the
pathways that include the perturbed reaction $i$. However, to precisely
conduct this summation, we need to know the whole stoichiometry matrix or
metabolic network. Hence, another relation closed only for the reaction fluxes
${{\bf v}}$ is required for application without the need to know the details
of the metabolic systems.

To derive such a relation, we consider effective changes in the stoichiometric
coefficients $S_{\alpha i}$ for reaction $i$ as metabolic inhibitions: e.g.,
inhibition of enzymes, administration of metabolite analogs, leakage of
metabolites, and inefficiency in the allocation of some resource. We then
obtain an equality on the optimized reaction fluxes $\hat{{{\bf v}}}$,
formally similar to Eq. (3) (see SM, Sec. S4 for derivation):

$\displaystyle\frac{\partial\hat{v}_{i}(S,{{\bf I}})}{\partial
q^{\nu}_{i}}=-\hat{v}_{i}(S,{{\bf I}})\frac{\partial\hat{v}_{i}(S,{{\bf
I}})}{\partial I_{\nu}}$ (4)

by defining the metabolic price $q^{\nu}_{i}$ of reaction $i$ in terms of
$\nu$ as a function of $S$, instead of the metabolic price
$p^{\nu}_{i^{\prime}}$ of pathway ${i^{\prime}}$ as a function of $K$,

$\displaystyle
q^{\nu}_{i}:=\sum_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}}-S_{\alpha
i}\left.\frac{\partial\hat{v}_{i}}{\partial
I_{\alpha}}\middle/\frac{\partial\hat{v}_{i}}{\partial I_{\nu}}\right..$ (5)

The coefficient ${\left.\frac{\partial\hat{v}_{i}}{\partial
I_{\alpha}}\middle/\frac{\partial\hat{v}_{i}}{\partial
I_{\nu}}\right.}\,{=:}\,c^{\nu}_{\alpha}(i)$ quantifies the number of units of
component $\nu$ that can compensate for one unit of $\alpha$ in reaction $i$
and is experimentally measurable. For example, if $\nu$ is glucose and
$\alpha$ is another metabolite such as an amino acid, $c^{\nu}_{\alpha}$
indicates how many units of glucose are required to compensate for one unit of
the amino acid, similar to the “glucose cost” in previous studies
chen2022yeast .

For the linear response relation (4), it is sufficient to calculate only the
change in metabolic price (not the metabolic price itself), which depends on
the type of manipulations in concern: (I) manipulations leading to the loss of
a single component and (II) those leading to the loss of multiple components.

If experimental manipulation causes the loss of a single component
$\alpha\,(\in{\mathcal{M}}\cup{\mathcal{C}})$ in reaction $i$, $S_{\alpha i}$
effectively changes only for that $\alpha$ [Fig. 3(a)]. In such a case, the
metabolic price change is just given by $\Delta q^{\alpha}_{i}\,{=}\,\Delta
S_{\alpha i}$. An example of such experimental manipulations is the
administration of an analog to a reactant of a multibody reaction: if $\alpha$
and $\beta$ react [see Fig. 3(a)], the metabolic analog of $\beta$ can produce
incorrect metabolite(s) with $\alpha$, leading to the loss of $\alpha$, and
thus, reaction $i$ requires more $\alpha$ to produce the same number of
products, causing effective increases in the input stoichiometric coefficient
$|S_{\alpha i}|$. Another example is the changes in the total capacity and
effective stoichiometry for a constraint: for example, the mitochondrial
volume capacity will work as such a constraint and can be genetically
manipulated malina2021adaptations ; van2001modulating ; raghevendran2006hap4 .
Equation (4) for case (I) is numerically confirmed in Fig. 3(a).

Figure 3: Responses of the optimized reaction fluxes $\hat{{{\bf v}}}$ against
metabolic inhibitions of reaction $i$. As the simplest example, a two-body
reaction of components $\alpha$ and $\beta$ is illustrated in the upper
panels. (a) Metabolic inhibition of a single component ($\alpha$ in upper
panel): case (I). The horizontal axis shows the responses to the available
amount of a constraint $\alpha\,(\in{\mathcal{C}})$,
$\hat{v}_{i}\Delta\hat{v}_{i}(\Delta I_{\alpha})/\Delta I_{\alpha}$, and the
vertical axis does those to metabolic inhibitions, $\Delta\hat{v}_{i}(\Delta
S_{\alpha i})/\Delta q^{\alpha}_{i}$. (b) Metabolic inhibition of multiple
components ($\alpha$ and $\beta$ in upper panel): case (II). Responses of the
reaction flux $\hat{v}_{i}$ to effective changes in the input stoichiometric
coefficients $\\{S_{\mu i}\\}_{\mu\in{\mathcal{M}}}$ for reaction $i$,
$\Delta\hat{v}_{i}(\\{\Delta S_{\mu i}\\}_{\mu\in{\mathcal{M}}})/\Delta
q^{\mathrm{Glc}}_{i}$, are plotted against those to intake changes,
$\hat{v}_{i}\Delta\hat{v}_{i}(\Delta I_{\mathrm{Glc}})/\Delta
I_{\mathrm{Glc}}$. $I_{\mathrm{Glc}}\,{=}\,8.05\,\mathrm{[mmol/gDW/h]}$. Each
marker denotes a different reaction $i$.

If metabolic inhibition of multiple reactant species of a reaction $i$ is
simultaneously caused, the stoichiometric coefficients $S_{\mu i}$ for
multiple reactants $\mu$ of reaction $i$ will effectively change [Fig. 3(b)].
Accordingly, the reaction price changes by $\Delta
q^{\nu}_{i}\,{=}\,\sum_{\mu}\Delta S_{\mu i}c^{\nu}_{\mu}(i)$. In experiments,
such cases would correspond to the inhibition of enzymes, leakage of the
intermediate complex of reaction $i$, and so forth. Even in this case (II),
the linear relation (4) is verified by numerically calculating the price
changes of reaction $j$ defined in Eq. (5) with the E.coli core model
including 77 reactions [Fig. 3(b)] as well as a larger-scale metabolic model
including 931 reactions reed2003expanded  (SM, Fig. S2). Here, although the
precise calculation of the coefficients $c^{\nu}_{\mu}(i)$ requires
information regarding not only the responses of $\hat{v}_{i}$ to $I_{\nu}$ but
also those to $I_{\mu}$, they can be approximated in ways easier and
independent of reaction $i$. For example, under extreme situations in which
only the carbon sources limit the objective reaction, $c^{\nu}_{\mu}$ should
be the ratios of the carbon numbers of species $\mu$ and $\nu$; alternatively,
the simplest approximation could be just taking $c^{\nu}_{\mu}$ as unities.
Even with these approximations, the relation (4) appears to hold well (SM,
Fig. S3), and thus such approximations will be useful for qualitatively
predicting whether metabolic inhibition promotes or suppresses the reaction of
interest.

Remarkably, our above argument does not depend on specific choices of
objective reaction $o$, whereas we have utilized the biomass synthesis
reaction as $o$ (Figs. 2 and 3). To highlight the independence of the relation
(4) from cellular objective functions, we also numerically confirmed that it
is satisfied even when objective reaction $o$ is set as a reaction for matter
production, such as ethanol or ATP synthesis (SM, Fig. S4). These synthesis
reactions are often considered as the objectives for metabolically engineered
cells stephanopoulos1998metabolic ; schuetz2007systematic ;
gianchandani2008predicting .

In the present study, we showed that the metabolic responses against resource
availability and those against metabolic inhibitions are negatively
proportional. The quantitative relations we found should be universally
satisfied with arbitrary reaction networks, constraints, and objective
functions of cells. In particular, although the predicted optimal metabolic
states can drastically depend on the assumed objective function, the relations
of the responses should be always satisfied independent of it (see also SM,
Fig. S4). Even though we can never know the true objective function of cells,
we can still predict the metabolic responses.

In the linear relations, the metabolic responses against different
perturbations are linked because both are determined from the identical
objective function and constraints (see also Fig. 1). It is similar to the
linear response theories in statistical mechanics: they are derived from the
fact that different thermodynamic quantities are given as the derivatives of
an identical thermodynamic potential LANDAU1980333 . Note here that the
thermodynamic potential works as an objective function: e.g., entropy is
maximized at the thermal equilibrium.

The independence from cellular objective functions is derived from the
microeconomic formulation for metabolic regulation and application of the
Slutsky equation in economics. Note that the Slutsky equation basically
requires detailed information regarding the objective functions (utility
functions in economics) because it includes a term for the so-called
substitution effect that quantifies the substitutability of goods and depends
on the objective functions (see also SM, Sec. S4). However, the term
disappears when applied to metabolism because the law of mass conservation
implies the non-substitutability of reactants.

Although the linear relations (3-4) are general due to the generality of the
law of mass conservation, there are also some limitations. First, since our
results rely on the assumption of optimal metabolic regulation, they will not
hold in suboptimal metabolic responses; conversely, any observed deviation
from the relations (3-4) will indicate the suboptimality in the regulation of
the real metabolic system in question. Second, our linear relations work only
for continuous metabolic responses. Third, the approximation of coefficients
$c_{\mu}^{\nu}$ in Eqs. (4-5) could be prohibitive when the coefficients
become negative, e.g., in the case the fluxes from different nutrient sources
must be balanced for a metabolic reaction of interest and an increase in one
source promotes the reaction while an increase in another source inhibits it.

Because our results are valid regardless of how abstract the concerned model
is, from coarse-grained toy models to genome-scale metabolic networks, they
would be important both for quantitative predictions and for discovering
qualitatively novel phenomena. The Warburg effect or overflow metabolism is a
prominent example of the latter. In the Warburg effect, as the amount of the
carbon source taken up by a cell increases, the cell decreases the flux of the
respiration pathway and utilize fermentation or aerobic glycolysis instead
Heiden-Warburg ; OMbook . From the relation (3), one can immediately predict
that the inhibition of respiration (e.g., administration of uncouplers of
respiration WeakAcid2 ) will counterintuitively increase the respiration flux.
Such an increase in the respiration flux was observed in a coarse-grained
model, which was termed the drug-induced reverse Warburg effect
yamagishi2021microeconomics . This phenomenon has been indeed reported in
several published experiments WeakAcid ; WeakAcid2 ; WeakAcid3 ; Therapies-
induced-reversed-Warburg .

Likewise, for controlling cellular metabolic states, e.g., for metabolic
engineering and medicine, some counterintuitive manipulations can promote
pathway or reaction fluxes. Although metabolic inefficiency is considered to
suppress the flux in general, when an increase in the intake of a substrate
suppresses a pathway or reaction flux, making the metabolic pathway or
reaction less efficient will counterintuitively promote the flux (see also SM,
Sec. S2 and Fig. S1 for an example of coarse-grained models). In experimental
application, the intake or total capacity can be altered by shifts in
environmental conditions, genetic manipulations, and so forth. Changes in the
metabolic price can be also implemented in various ways: e.g., administration
of a metabolite analog, leakage of a metabolite, addition of an alternative
pathway or reaction through metabolic engineering manipulations, and
inhibition of some enzyme that will lead to the accumulation of the reactants
and possibly promote their excretion or conversion to other chemicals. They
cause a loss of reactants, and thus, the corresponding reaction(s) require
more metabolites to produce the same number of products.

The relations (3-4) allow us to predict the responses of an arbitrary reaction
or pathway flux to metabolic inhibitions only by measuring its fluxes in
several different nutrient conditions, and vice versa. The predictions do not
require detailed information regarding the concerned intracellular reaction
networks, and they are valid even when the precise estimation of effective
changes in the stoichiometric coefficients is difficult, at least
qualitatively (SM, Fig. S3). Therefore, they will be useful as quantitative
and qualitative guidelines to operate the metabolic states toward the
desirable directions in various fields such as microbiology, metabolic
engineering, and medicine.

The associated Python code and data to reproduce figures in this work are
available 111https://github.com/JFYamagishi/yamagishi-hatakeyama-2023.

We would like to acknowledge Wolfram Liebermeister, Chikara Furusawa, Yasushi
Okada, Kunihiko Kaneko, and Takuma Ōnishi for helpful discussions and useful
comments. This work was partially supported by JSPS KAKENHI Grant Numbers
JP21J22920 and JP21K15048.

## References

  * (1) E. Ilker, M. Hinczewski, Modeling the growth of organisms validates a general relation between metabolic costs and natural selection, Physical Review Letters 122 (23) (2019) 238101.
  * (2) S. Kar, D. S. Ray, Collapse and revival of glycolytic oscillation, Physical Review Letters 90 (23) (2003) 238102.
  * (3) J. F. Yamagishi, N. Saito, K. Kaneko, Advantage of leakage of essential metabolites for cells, Physical Review Letters 124 (4) (2020) 048101.
  * (4) D. Zwicker, R. Seyboldt, C. A. Weber, A. A. Hyman, F. Jülicher, Growth and division of active droplets provides a model for protocells, Nature Physics 13 (4) (2017) 408–413.
  * (5) G. Stephanopoulos, A. A. Aristidou, J. Nielsen, Metabolic engineering: principles and methodologies, Elsevier, Amsterdam, 1998.
  * (6) P. Murima, J. D. McKinney, K. Pethe, Targeting bacterial central metabolism for drug development, Chemistry & Biology 21 (11) (2014) 1423–1432.
  * (7) M. G. Vander Heiden, Targeting cancer metabolism: a therapeutic window opens, Nature Reviews Drug Discovery 10 (9) (2011) 671–684.
  * (8) U. E. Martinez-Outschoorn, M. Peiris-Pagés, R. G. Pestell, F. Sotgia, M. P. Lisanti, Cancer metabolism: a therapeutic perspective, Nature Reviews Clinical Oncology 14 (1) (2017) 11–31.
  * (9) B. Palsson, Systems Biology, Cambridge university press, 2015.
  * (10) E. Klipp, W. Liebermeister, C. Wierling, A. Kowald, Systems Biology: a textbook, John Wiley & Sons, New Jersey, 2016.
  * (11) P. B. Warren, J. L. Jones, Duality, thermodynamics, and the linear programming problem in constraint-based models of metabolism, Physical Review Letters 99 (10) (2007) 108101.
  * (12) R. Heinrich, S. Schuster, The modelling of metabolic systems. structure, control and optimality, Biosystems 47 (1-2) (1998) 61–77.
  * (13) A. Vazquez, Overflow metabolism: from yeast to marathon runners, Academic Press, London, 2017.
  * (14) V. A. Portnoy, D. Bezdan, K. Zengler, Adaptive laboratory evolution—harnessing the power of biology for metabolic engineering, Current Opinion in Biotechnology 22 (4) (2011) 590–594.
  * (15) H. P. Bonarius, G. Schmid, J. Tramper, Flux analysis of underdetermined metabolic networks: the quest for the missing constraints, Trends in Biotechnology 15 (8) (1997) 308–314.
  * (16) K. Raman, N. Chandra, Flux balance analysis of biological systems: applications and challenges, Briefings in Bioinformatics 10 (4) (2009) 435–449.
  * (17) B. Schnitzer, L. Österberg, M. Cvijovic, The choice of the objective function in flux balance analysis is crucial for predicting replicative lifespans in yeast, Plos One 17 (10) (2022) e0276112.
  * (18) B. Niebel, S. Leupold, M. Heinemann, An upper limit on gibbs energy dissipation governs cellular metabolism, Nature Metabolism 1 (1) (2019) 125–132.
  * (19) L. Onsager, Reciprocal relations in irreversible processes. i., Physical review 37 (4) (1931) 405.
  * (20) R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, Journal of the Physical Society of Japan 12 (6) (1957) 570–586.
  * (21) M. S. Green, Markoff random processes and the statistical mechanics of time-dependent phenomena. ii. irreversible processes in fluids, The Journal of chemical physics 22 (3) (1954) 398–413.
  * (22) H. R. Varian, Microeconomic analysis, WW Norton, New York, 1992.
  * (23) K. J. Lancaster, A new approach to consumer theory, Journal of political economy 74 (2) (1966) 132–157.
  * (24) J. F. Yamagishi, T. S. Hatakeyama, Microeconomics of metabolism: The warburg effect as giffen behaviour, Bulletin of Mathematical Biology 83 (12) (2021) 1–17.
  * (25) M. Scott, T. Hwa, Bacterial growth laws and their applications, Current Opinion in Biotechnology 22 (4) (2011) 559–565.
  * (26) M. Basan, S. Hui, H. Okano, Z. Zhang, Y. Shen, J. R. Williamson, T. Hwa, Overflow metabolism in escherichia coli results from efficient proteome allocation, Nature 528 (7580) (2015) 99–104.
  * (27) A. Vazquez, J. Liu, Y. Zhou, Z. N. Oltvai, Catabolic efficiency of aerobic glycolysis: the warburg effect revisited, BMC Systems Biology 4 (1) (2010) 1–9.
  * (28) M. Szenk, K. A. Dill, A. M. de Graff, Why do fast-growing bacteria enter overflow metabolism? testing the membrane real estate hypothesis, Cell Systems 5 (2) (2017) 95–104.
  * (29) C. Liao, T. Wang, S. Maslov, J. B. Xavier, Modeling microbial cross-feeding at intermediate scale portrays community dynamics and species coexistence, PLoS Computational Biology 16 (8) (2020) e1008135.
  * (30) A. Roy, D. Goberman, R. Pugatch, A unifying autocatalytic network-based framework for bacterial growth laws, Proceedings of the National Academy of Sciences 118 (33) (2021) e2107829118.
  * (31) M. Terzer, J. Stelling, Large-scale computation of elementary flux modes with bit pattern trees, Bioinformatics 24 (19) (2008) 2229–2235.
  * (32) J. D. Orth, R. M. Fleming, B. Ø. Palsson, Reconstruction and use of microbial metabolic networks: the core escherichia coli metabolic model as an educational guide, EcoSal plus 4 (1) (2010).
  * (33) C. H. Schilling, D. Letscher, B. Ø. Palsson, Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective, Journal of Theoretical Biology 203 (3) (2000) 229–248.
  * (34) Y. Chen, J. Nielsen, Yeast has evolved to minimize protein resource cost for synthesizing amino acids, Proceedings of the National Academy of Sciences 119 (4) (2022).
  * (35) C. Malina, R. Yu, J. Björkeroth, E. J. Kerkhoven, J. Nielsen, Adaptations in metabolism and protein translation give rise to the crabtree effect in yeast, Proceedings of the National Academy of Sciences 118 (51) (2021).
  * (36) A. J. van Maris, B. M. Bakker, M. Brandt, A. Boorsma, M. J. Teixeira de Mattos, L. A. Grivell, J. T. Pronk, J. Blom, Modulating the distribution of fluxes among respiration and fermentation by overexpression of hap4 in saccharomyces cerevisiae, FEMS Yeast Research 1 (2) (2001) 139–149.
  * (37) V. Raghevendran, K. R. Patil, L. Olsson, J. Nielsen, Hap4 is not essential for activation of respiration at low specific growth rates in saccharomyces cerevisiae, Journal of Biological Chemistry 281 (18) (2006) 12308–12314.
  * (38) J. L. Reed, T. D. Vo, C. H. Schilling, B. O. Palsson, An expanded genome-scale model of escherichia coli k-12 (ijr904 gsm/gpr), Genome Biology 4 (9) (2003) 1–12.
  * (39) R. Schuetz, L. Kuepfer, U. Sauer, Systematic evaluation of objective functions for predicting intracellular fluxes in escherichia coli, Molecular Systems Biology 3 (1) (2007) 119.
  * (40) E. P. Gianchandani, M. A. Oberhardt, A. P. Burgard, C. D. Maranas, J. A. Papin, Predicting biological system objectives de novo from internal state measurements, BMC Bioinformatics 9 (1) (2008) 1–13.
  * (41) L. Landau, E. Lifshitz, Chapter xii - fluctuations, in: Statistical Physics, 3rd Edition, Butterworth-Heinemann, Oxford, 1980, Ch. 12, pp. 333–400.
  * (42) M. G. Vander Heiden, L. C. Cantley, C. B. Thompson, Understanding the warburg effect: the metabolic requirements of cell proliferation, Science 324 (5930) (2009) 1029–1033.
  * (43) C. Verduyn, E. Postma, W. A. Scheffers, J. P. Van Dijken, Effect of benzoic acid on metabolic fluxes in yeasts: a continuous-culture study on the regulation of respiration and alcoholic fermentation, Yeast 8 (7) (1992) 501–517.
  * (44) E. Postma, C. Verduyn, A. W. Scheffers, J. P. Van Dijken, Enzymic analysis of the crabtree effect in glucose-limited chemostat cultures of saccharomyces cerevisiae., Applied and Environmental Microbiology 55 (2) (1989) 468–477.
  * (45) M. Gallmetzer, W. Burgstaller, Efflux of organic acids in penicillium simplicissimum is an energy-spilling process, adjusting the catabolic carbon flow to the nutrient supply and the activity of catabolic pathways, Microbiology 148 (4) (2002) 1143–1149.
  * (46) J. da Veiga Moreira, M. Hamraz, M. Abolhassani, L. Schwartz, M. Jolicœur, S. Pérès, Metabolic therapies inhibit tumor growth in vivo and in silico, Scientific reports 9 (1) (2019) 1–10.
  * (47) https://github.com/JFYamagishi/yamagishi-hatakeyama-2023.
  * (48) T. J. Clement, E. B. Baalhuis, B. Teusink, F. J. Bruggeman, R. Planqué, D. H. de Groot, Unlocking elementary conversion modes: Ecmtool unveils all capabilities of metabolic networks, Patterns 2 (1) (2021) 100177.
  * (49) A. Flamholz, E. Noor, A. Bar-Even, W. Liebermeister, R. Milo, Glycolytic strategy as a tradeoff between energy yield and protein cost, Proceedings of the National Academy of Sciences 110 (24) (2013) 10039–10044.
  * (50) E. Reznik, P. Mehta, D. Segrè, Flux imbalance analysis and the sensitivity of cellular growth to changes in metabolite pools, PLoS Computational Biology 9 (8) (2013) e1003195.

Supplemental Material

##  S1 Equivalence between constraint-based modeling (CBM) and microeconomic
formulations

CBM formulation: optimization problems with reaction fluxes ${{\bf v}}$ as
variables. In the framework of CBM in systems biology, intracellular metabolic
regulation is formulated as linear programming (LP) problems in which the
variables are the fluxes ${{\bf v}}$ of reactions. As discussed below, LP
problems in CBM are generally equivalent to optimization problems in the
microeconomic theory of consumer choice.

By breaking down each reversible reaction into two irreversible reactions
(i.e., its forward and backward components), a non-negative
$|\mathcal{R}|$-dimensional vector ${{\bf
v}}:=\\{v_{i}\\}_{i\in{\mathcal{R}}}$ represents the fluxes of all reactions.
Then, by assuming the stationarity of the intracellular concentrations of non-
exchangeable species $\mu\;(\in{\mathcal{M}}\backslash\mathcal{E})$, a general
formulation of CBM palsson2015systems ; klipp2016systems  is given as follows:

$\displaystyle\max_{{{\bf v}}\geq{\bf 0}}$
$\displaystyle\;v_{o}\;\;\mathrm{s.t.}$ (S2)
$\displaystyle\sum_{j\in{\mathcal{R}}}S_{\mu
j}v_{j}=0\;\;(\mu\in{\mathcal{M}}\backslash\mathcal{E})$
$\displaystyle\sum_{j\in{\mathcal{R}}}S_{\alpha j}v_{j}+I_{\alpha}\geq
0\;\;(\alpha\in\mathcal{E}\cup{\mathcal{C}})$

Because $\sum_{j\in{\mathcal{R}}}S_{\mu j}v_{j}$ is equal to the excess
production of species $\mu$, Eq. (S2) represents that the production and
degradation of internal metabolites must be balanced. With respect to Eq.
(S2), if $\alpha$ is a species ($\alpha\in{\mathcal{E}}\subset{\mathcal{M}}$),
the corresponding inequality, $\sum_{j\in{\mathcal{R}}}S_{\alpha
j}v_{j}+I_{\alpha}\geq 0$, represents that exchangeable species $\alpha$ with
intake $I_{\alpha}>0$ are taken in and species $\alpha$ with efflux
$I_{\alpha}<0$ are (forcibly) leaked or degraded; in contrast, if $\alpha$ is
a constraint ($\alpha\in{\mathcal{C}}$), it represents a non-stoichiometric
constraint, e.g., allocation of some limited resource such as proteins
scott2011bacterial ; OM , intracellular space vazquez2010catabolic ; OMbook ,
membrane surface memRealEstate2 , and Gibbs energy dissipation Heinemann-
Gibbs2019 . Stoichiometric coefficient $S_{\alpha i}$ for some constraint
$\alpha\,(\in{\mathcal{C}})$ is typically negative or zero but can also be
positive. Note that ${\mathcal{C}}$ can also include other constraints like
the upper and lower bounds of the flux of reaction $i$.

Optimization problems with pathway fluxes ${{\bf f}}$ as variables. The usual
CBM formulation (S2-S2) with reaction fluxes ${{\bf v}}$ as variables is
equivalent to another LP problem with pathway fluxes ${{\bf f}}$ as variables
(see also Fig. 1 in the main text).

Here, a metabolic pathway is a linked sequence of reactions. Pathway fluxes
${{\bf f}}$ are related to reaction fluxes as ${{\bf v}}=P{{\bf f}}$, with a
pathway matrix $P:=\\{P_{i{i^{\prime}}}\geq
0\;|\;i\in{\mathcal{R}},{i^{\prime}}\in{\mathcal{P}}\\}$ which represents that
pathway $iprime$ comprises $P_{i{i^{\prime}}}$ units of reaction $i$. The
stoichiometry matrix $K$ for pathways expresses the metabolic pathways in
“species space” and is related to the stoichiometry matrix $S$ for reactions
with $K=SP$.

When $P$ is taken as (linear combinations of) elementary flux modes (EFMs;
i.e., extreme rays of the flux cone), Eq. (S2) is autonomously satisfied (and
matrix $K:=SP$ is then related to so-called elementary conversion modes
clement2021unlocking ). Then, the LP problem (S2-S2) for CBM can be rewritten
into another LP problem with pathway fluxes ${{\bf f}}$ as variables:

$\displaystyle\max_{\Lambda,{{\bf f}}\geq{\bf
0}}\;\Lambda\quad\mathrm{s.t.}\quad
I_{\alpha}+\sum_{j^{\prime}\in{\mathcal{P}}}K_{\alpha
j^{\prime}}f_{j^{\prime}}\geq-S_{\alpha
o}\Lambda\;\;(\alpha\in\mathcal{E}\cup{\mathcal{C}})$ (S3)

where the components required for objective reaction $o\,(\in{\mathcal{R}})$
are termed as objective components
${\mathcal{O}}:=\\{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}\;|\;S_{\alpha
o}<0\\}$.

Derivation of the microeconomic formulation. The min function for the
objective function, $\Lambda:=\min(A,B)$, can be replaced by constraints of
$\Lambda\leq A$ and $\Lambda\leq B$, and vice versa. Therefore, the
optimization problem (S3) is equivalent to the optimization of the Leontief
utility function [Eq. (1) in the main text]—the minimum of multiple
“complementary” objectives—under the constraints of Eq. (2) in the main text.
In contrast, the microeconomic formulation of optimization problems (Eqs.
(1-2) in the main text) can be converted to LP problems in the form of Eq.
(S3), and thus, they are equivalent.

##  S2 A simple example of microeconomic formulation of metabolic regulation

As a concrete example, we here introduce a simple, analytically-solvable
metabolic model (see also Fig. S1). It can be the simplest model of a
metabolic system in which multiple pathways can produce a common metabolite
from a single resource, e.g., the co-utilization of respiration and
fermentation or that of Embden–Meyerhoff–Parnass (EMP) and Entner–Doudoroff
(ED) glycolytic pathways EMP/ED .

This model consists of $2$ metabolites, ${\mathcal{M}}=\\{m_{1},m_{2}\\}$, and
$1$ constraint ${\mathcal{C}}=\\{\rho\\}$, as well as $3$ reactions,
${\mathcal{R}}=\\{A,B$, objective$\\}$. Note that $\rho$ can be any kind of
limited resource, e.g., such as the intracellular volume or solvent capacity
vazquez2010catabolic ; OMbook , the total amount of proteins OM , and the
total area of membrane surfaces memRealEstate2 . We assume that $m_{1}$ is,
out of ${\mathcal{M}}$, the only exchangeable species, i.e.,
${\mathcal{E}}=\\{m_{1}\\}$, reactions $A,B$ both produce $m_{2}$ by consuming
$m_{1}$ and $\rho$, and the objective reaction $o$ requires $m_{2}$ and
$\rho$. The stoichiometry matrix,
$S\in\mathbb{R}^{({\mathcal{M}}\cup{\mathcal{C}})\times{\mathcal{R}}}$, is
then given as

$\displaystyle\begin{split}S=\begin{pmatrix}-p^{m_{1}}_{A}&-p^{m_{1}}_{B}&0\\\
S_{m_{2}A}&S_{m_{2}B}&-1\\\
-p^{\rho}_{A}&-p^{\rho}_{B}&-1\end{pmatrix},\end{split}$

with positive parameters (i.e., “price” of metabolic pathways)
$p^{\alpha}_{i}$ $(\alpha=m_{1},\rho;\;i=A,B)$. For simplicity, the
stoichiometric coefficients are normalized as $|S_{m_{2}o}|=|S_{\rho o}|=1$,
whereas it does not lose its generality. With this matrix, the LP problem for
CBM with reaction fluxes ${{\bf v}}$ as variables is represented as

$\displaystyle\max_{{{\bf
v}}\geq\textbf{0}}\;v_{\mathrm{BM}}\quad\mathrm{s.t.}\quad$
$\displaystyle\sum_{j}S_{\mu j}v_{j}+I_{\mu}\geq 0\quad(\mu=m_{1},\rho),$
$\displaystyle\sum_{j}S_{m_{2},j}v_{j}=0.$

In this simple model, reactions $A,B$ are metabolic pathways as they are,
since they connect the nutritional species $m_{1}$ to objective component
$m_{2}$. Then, we do not distinguish the fluxes of reactions ${{\bf v}}$ and
those of metabolic pathways ${{\bf f}}$. Thus, the microeconomic formulation
of the above LP problem is the maximization of

$\displaystyle\Lambda({{\bf
f}}):=\min\left[S_{m_{2}A}f_{A}+S_{m_{2}B}f_{B},\;I_{\rho}-p^{\rho}_{A}f_{A}-p^{\rho}_{B}f_{B}\right]\quad\mathrm{s.t.}\quad
p^{m_{1}}_{A}f_{A}+p^{m_{1}}_{B}f_{B}\leq I_{m_{1}},$

and the relations (3) and (4) in the main text are equivalent in this model:

$\displaystyle\frac{\partial\hat{f}_{i}}{\partial
p^{m_{1}}_{j}}=-\hat{f}_{j}\frac{\partial\hat{f}_{i}}{\partial
I_{m_{1}}}.\quad(i,j=A,B)$ (S4)

We can analytically obtain the optimized solutions for this model, from which
we can directly confirm that the relation (S4) is indeed satisfied as follows.
Since the solutions depend on the relative values of the parameters, we
analyze the behavior of the model here by dividing it into two cases: whether
some reaction flux is suppressed by the increase in nutrient intake
$I_{m_{1}}$ or not.

When no reaction is suppressed by increasing the nutrient intake. Without any
trade-off between reactions $A,B$, only the more efficient reaction is used in
the optimized solutions. In such case, no reaction flux is suppressed by
increasing the nutrient intake and metabolic inhibition of a reaction always
causes its suppression.

For example, when reaction $A$ is more efficient, the optimal flux $\hat{{{\bf
f}}}$ is just given by

$\displaystyle\left(\hat{f}_{A},\hat{f}_{B}\right)=\left(\min\left[\frac{I_{m_{1}}}{p^{m_{1}}_{A}},\frac{I_{\rho}}{p^{\rho}_{A}}\right],0\right),$
(S5)

and the flux $\hat{f}_{A}$ is promoted or sustained by increasing $I_{m_{1}}$.
Also, from Eq. (S5), when $I_{m_{1}}/p^{m_{1}}_{A}\leq I_{\rho}/p^{\rho}_{A}$
holds, the metabolic responses are given as

$\displaystyle\frac{\partial\hat{f}_{A}}{\partial
I_{m_{1}}}=\frac{1}{p^{m_{1}}_{A}},\quad\frac{\partial\hat{f}_{A}}{\partial
p^{m_{1}}_{A}}=-\frac{I_{m_{1}}}{(p^{m_{1}}_{A})^{2}},\quad\frac{\partial\hat{f}_{A}}{\partial
p^{m_{1}}_{A}}=-\hat{f}_{A}\frac{\partial\hat{f}_{A}}{\partial I_{m_{1}}},$

and

$\displaystyle\hat{f}_{A}\frac{\partial\hat{f}_{B}}{\partial
I_{m_{1}}}=\frac{\partial\hat{f}_{B}}{\partial
p^{m_{1}}_{A}}=0,\quad\hat{f}_{B}\frac{\partial\hat{f}_{i}}{\partial
I_{m_{1}}}=\frac{\partial\hat{f}_{i}}{\partial p^{m_{1}}_{B}}=0\;\;(i=A,B),$

and thus the relation (S4) holds. In the case
$I_{m_{1}}/p^{m_{1}}_{A}>I_{\rho}/p^{\rho}_{A}$, the optimal fluxes are
independent of the value of $I_{m_{1}}$ and the relation (S4) is evidently
satisfied as

$\displaystyle\hat{f}_{j}\frac{\partial\hat{f}_{i}}{\partial
I_{m_{1}}}=\frac{\partial\hat{f}_{i}}{\partial p^{m_{1}}_{j}}=0\;\;(i,j=A,B).$

When reaction $A$ is suppressed by increasing the nutrient intake. With some
trade-off between reactions $A,B$, the situation is non-trivial: the co-
utilization of different pathways $A,B$ and switching between them can be
observed. In particular, we here assume a trade-off $s_{m_{2}A}>s_{m_{2}B}$
and $p^{\rho}_{A}>p^{\rho}_{B}$, i.e., reaction $A$ is more efficient in
producing $m_{2}$ but requires more amount of constraint $\rho$ than the
alternative reaction $B$. When the trade-off matters, the optimized flux of
reaction $A$ is suppressed by increasing nutrient intake $I_{m_{1}}$ and it is
promoted by its metabolic inefficiency.

Figure S1:  A simple example of microeconomic formulation of metabolic
regulation.

We also assume that $p^{m_{1}}_{A}$ is not so large compared to
$p^{m_{1}}_{B}$ that $p^{m_{1}}_{B}B_{0}>p^{m_{1}}_{A}A_{0}$ holds, where
$A_{0}:=I_{\rho}/(S_{m_{2}A}+p^{\rho}_{A})$ and
$B_{0}:=I_{\rho}/(S_{m_{2}B}+p^{\rho}_{B})$. The optimal fluxes $\hat{{{\bf
f}}}$ is then calculated as

$\displaystyle\left(\hat{f}_{A},\hat{f}_{B}\right)=\begin{cases}\left(\frac{I_{m_{1}}}{p^{m_{1}}_{A}},0\right)\quad\mathrm{if}\;I_{m_{1}}\leq
p^{m_{1}}_{A}A_{0}\\\
\left(\frac{p^{m_{1}}_{B}B_{0}-I_{m_{1}}}{p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0}}A_{0},\frac{I_{m_{1}}-p^{m_{1}}_{A}A_{0}}{p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0}}B_{0}\right)\quad\mathrm{if}\;p^{m_{1}}_{B}B_{0}\geq
I_{m_{1}}\geq p^{m_{1}}_{A}A_{0}\\\
\left(0,B_{0}\right)\quad\mathrm{if}\;I_{m_{1}}\geq
p^{m_{1}}_{B}B_{0}\end{cases}$

From the dependence of this solution on $I_{m_{1}}$ and $p^{m_{1}}_{i}$, we
can immediately confirm that the relation (S4) is satisfied.

In the case of $I_{m_{1}}\leq p^{m_{1}}_{A}A_{0}$ or $I_{m_{1}}\geq
p^{m_{1}}_{B}B_{0}$, either reaction is used and thus the relation (S4) is
evidently satisfied as in the above case.

In contrast, with the medium value of the intake $I_{m_{1}}$ of metabolite
$m_{1}$ that satisfies $p^{m_{1}}_{B}B_{0}\geq I_{m_{1}}\geq
p^{m_{1}}_{A}A_{0}$, the metabolic responses against $I_{m_{1}}$ are:

$\displaystyle\frac{\partial\hat{f}_{A}}{\partial I_{m_{1}}}$ $\displaystyle=$
$\displaystyle-\frac{A_{0}}{p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0}},\quad\frac{\partial\hat{f}_{B}}{\partial
I_{m_{1}}}=\frac{B_{0}}{p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0}}.$

Note here that, as $p^{m_{1}}_{B}B_{0}>p^{m_{1}}_{A}A_{0}$ holds from the
assumption,

$\displaystyle\frac{\partial\hat{f}_{A}}{\partial
I_{m_{1}}}<0,\;\frac{\partial\hat{f}_{B}}{\partial I_{m_{1}}}>0.$

Thus, we can immediately obtain the qualitative prediction that the inhibition
of reaction/pathway $A$ promotes its own flux $\hat{f}_{A}$ as well as the
quantitative prediction (S4). Since the price responses are

$\displaystyle\frac{\partial\hat{f}_{A}}{\partial
p^{m_{1}}_{A}}=\frac{p^{m_{1}}_{B}B_{0}-I_{m_{1}}}{(p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0})^{2}}(A_{0})^{2}=-\hat{f}_{A}\frac{\partial\hat{f}_{A}}{\partial
I_{m_{1}}},$ $\displaystyle\frac{\partial\hat{f}_{A}}{\partial
p^{m_{1}}_{B}}=\frac{I_{m_{1}}-p^{m_{1}}_{A}A_{0}}{(p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0})^{2}}A_{0}B_{0}=-\hat{f}_{B}\frac{\partial\hat{f}_{A}}{\partial
I_{m_{1}}},$ $\displaystyle\frac{\partial\hat{f}_{B}}{\partial
p^{m_{1}}_{A}}=-\frac{p^{m_{1}}_{B}B_{0}-I_{m_{1}}}{(p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0})^{2}}A_{0}B_{0}=-\hat{f}_{A}\frac{\partial\hat{f}_{B}}{\partial
I_{m_{1}}},$ $\displaystyle\frac{\partial\hat{f}_{B}}{\partial
p^{m_{1}}_{B}}=-\frac{I_{m_{1}}-p^{m_{1}}_{A}A_{0}}{(p^{m_{1}}_{B}B_{0}-p^{m_{1}}_{A}A_{0})^{2}}(B_{0})^{2}=-\hat{f}_{B}\frac{\partial\hat{f}_{B}}{\partial
I_{m_{1}}},$

the linear relation (S4) is indeed satisfied.

From the biological perspective, by regarding reactions/pathways $A$ and $B$
as respiration and fermentation pathways, respectively, the above model can be
interpreted as a minimal model for overflow metabolism or the Warburg effect
OMbook . In one interpretation, $m_{1}$, $m_{2}$, $\rho$, and $o$ can be
regarded as glucose, ATP, the solvent capacity, and biomass synthesis
reaction, respectively. Then, the trade-off and the suppression of respiration
against the increase in carbon intake are consistent with the empirical
observations. Using our linear relation, we can quantitatively predict the
stimulation of the respiration flux (i.e., the drug-induced reverse Warburg
effect) by the administration of uncouplers of respiration.

##  S3 Details of numerical experiments

In our numerical simulations, stoichiometric coefficients $S_{\alpha i}$ for
species $\alpha\in{\mathcal{M}}$ are given by the E.coli core model
orth2010reconstruction . $S_{\alpha i}$’s for constraints
$\alpha\in{\mathcal{C}}$ are randomly chosen from an interval $[-1,0]$ to show
that the results do not depend on such details; here, as a natural assumption,
we set $S_{\alpha i}$ to be identical for the pair of reactions split from the
same reversible reaction.

In most numerical simulations, the objective reaction $o$ is chosen as the
biomass synthesis reaction, whereas our results do not depend on the choice of
the objective reactions: e.g., synthesis reaction of ethanol or ATP (see also
Fig. S4).

We plotted only reactions with continuous price responses in Figs. 2 and 3 in
the main text and Figs. S2-S4. In our numerical calculations, the response to
price change is considered to be discontinuous if decreasing $\Delta p$
changes the reaction flux discontinuously at $\Delta p\to 0$.

On the uniqueness of the solutions. For simplicity, the arguments in the main
text are based on the premise of the uniqueness of the solution. This
assumption indeed holds true with suitable penalty terms explained below.
Furthermore, it must be biologically natural because we considered the
intracellular metabolic responses around a steady state here.

An LP can sometimes have multiple solutions. If an LP has multiple solutions
because of the existence of ineffective inequalities, it is natural to define
the unique solution as the solution with the least consumption rates of
species: i.e., slightly modify the objective function as
$v_{\mathrm{BM}}-\epsilon\sum_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}}\left(\sum_{i\in{\mathcal{R}}}-S_{\alpha
i}v_{i}\right)$ with small $\epsilon>0$ (here, note that
$\sum_{i\in{\mathcal{R}}}-S_{\alpha i}v_{i}$ is equal to the consumption rate
of $\alpha$, while $I_{\alpha}$ is its upper bound). In addition, for the pair
of two irreversible reactions formed by separating each reversible reaction,
the smaller of the two must be subtracted and set to zero for the uniqueness
of the solution by adding a penalty $-\epsilon\sum_{i}v_{i}$. Throughout the
numerical calculations in the main text, $\epsilon$ is basically set to
$\lesssim 10^{-5}$, whereas the results do not depend on $\epsilon$ as long as
$\epsilon$ is sufficiently small.

Leakage of intermediate species. To calculate the reaction price $q^{\nu}_{i}$
(Eq. (5) in the main text), we must calculate the responses of $\hat{v}_{i}$
against changes in $I_{\mu}$ for each reactant $\mu$ of reaction $i$. For
calculating the responses to $I_{\mu}$ of each $\mu$, we assumed that species
$\mu$ is exchangeable. Such changes do not alter the optimized solution
$\hat{{{\bf v}}}$ in the numerical calculations we conducted.

Choice of metabolic pathways in Fig. 2 in the main text. In this study, we
consider (linear sums of) EFMs or extreme pathways for the stoichiometry
without objective reaction $o$ as the metabolic pathways. That is, we consider
metabolic pathways from components $\alpha$ with influxes $I_{\alpha}>0$ to
objective components ${\mathcal{O}}$.

##  S4 Derivation of Eqs. (3-5) in the main text

Derivation of Eq. (3): Slutsky equation in microeconomics. When given a
utility function $u(\textbf{x})$, we define $\hat{x}_{i}({\bf p},I)$ as the
optimal demand for good $i$ determined as a function of the price of the good
${\bf p}$ and income $I$. By defining $E({\bf p},u)$ as “the minimum income
required to achieve a certain utility value $u$,” we can represent “the
minimum demand for good $i$ required to achieve a utility value $u$” as
$h_{i}({\bf p},u):=\hat{x}_{i}({\bf p},E({\bf p},u))$.

Differentiating this function $h_{i}({\bf p},u)$ with respect to $p_{j}$
yields

$\displaystyle\frac{\partial h_{i}({\bf p},\hat{u}({\bf p},I))}{\partial
p_{j}}=\frac{\partial\hat{x}_{i}({\bf p},I)}{\partial
p_{j}}+\frac{\partial\hat{x}_{i}({\bf p},I)}{\partial I}\frac{\partial E({\bf
p},\hat{u}({\bf p},I))}{\partial p_{j}},$

where $\hat{u}({\bf p},I)$ represents the maximum utility under a given price
${\bf p}$ and income $I$. Due to optimality, the last term ${\partial E({\bf
p},\hat{u}({\bf p},I))}/{\partial p_{j}}$ equals to the optimal demand
$\hat{x}_{j}({\bf p},I)$.

Accordingly, we obtain the Slutsky equation that describes the response of
demand $\hat{x}_{i}({\bf p},I)$ to changes in price $p_{j}$:

$\displaystyle\frac{\partial\hat{x}_{i}({\bf p},I)}{\partial
p_{j}}=\frac{\partial h_{i}({\bf p},\hat{u}({\bf p},I))}{\partial
p_{j}}-\hat{x}_{j}({\bf p},I)\frac{\partial\hat{x}_{i}({\bf p},I)}{\partial
I}.$ (S6)

The first term ${\partial h_{i}({\bf p},\hat{u})}/{\partial p_{j}}$ represents
the substitution effect caused by relative changes in the price of each good;
particularly, the “self-substitution effect” for $i=j$ is always non-positive
Varian . In contrast, the second term $\hat{x}_{j}({\bf
p},I)\frac{\partial\hat{x}_{i}({\bf p},I)}{\partial I}$ reflects the income
effect, which can be either positive or negative. This represents the effect
that an increase in the price of a good leads to an effective decrease in
income, which also changes the demand for goods.

The law of mass conservation in metabolism corresponds has the perfect
complementarity in economics. Accordingly, the substitution effect is always
null (at a kink) Lancaster1966 ; yamagishi2021microeconomics . Then, by noting
that $|\mathcal{P}|$ metabolic pathways correspond to goods in economics in
our mapping between metabolism and microeconomics, we immediately obtain Eq.
(3) in the main text from the Slutsky equation (S6).

Derivation of Eqs. (4-5) from Eq. (3) in the main text. Noting $\hat{{{\bf
v}}}=P\hat{{{\bf f}}}$ and multiplying Eq. (3) in the main text by pathway
matrix $P$ from the left, we immediately obtain equality for the responses of
reaction fluxes ${{\bf v}}$:

$\displaystyle\frac{\partial\hat{v}_{i}(K,{\bf I})}{\partial
p^{\nu}_{j^{\prime}}}=-f_{j^{\prime}}(K,{\bf
I})\frac{\partial\hat{v}_{i}(K,{\bf I})}{\partial
I_{\nu}}\quad(i\in{\mathcal{R}},\;j^{\prime}\in{\mathcal{P}}),$ (S7)

Assuming that a stoichiometric coefficient $S_{\mu j}$ of reaction $j$ is
effectively altered to $S_{\mu j}-\Delta S_{\mu j}$, the ${j^{\prime}}$-th
column vector of matrix $K:=SP$ changes as

$\displaystyle\\{K_{\alpha{j^{\prime}}}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}}\;$
$\displaystyle\to$ $\displaystyle\;\\{K_{\alpha{j^{\prime}}}-\Delta
K_{\alpha{j^{\prime}}}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}}:=\\{\sum_{k\in{\mathcal{R}}}S_{\alpha
k}P_{k{j^{\prime}}}-\delta_{\alpha\mu}\Delta S_{\mu
j}P_{j{j^{\prime}}}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}},$

with Kronecker delta $\delta_{\alpha\mu}$. That is, the metabolic price of
pathway ${j^{\prime}}$ for metabolite $\mu$ changes by $\Delta
p^{\mu}_{j^{\prime}}(\\{\delta_{\alpha\mu}\Delta S_{\mu
j}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}})=\Delta S_{\mu
j}P_{j{j^{\prime}}}$.

Therefore, if multiple stoichiometric coefficients are simultaneously altered,

$\displaystyle\Delta\hat{v}_{i}(\\{\Delta S_{\alpha
j}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}})$ $\displaystyle=$
$\displaystyle\sum_{\mu\in{\mathcal{M}}}\sum_{{j^{\prime}}\in{\mathcal{P}}}\frac{\partial\hat{v}_{i}(K,{\bf
I})}{\partial p^{\mu}_{j^{\prime}}}\Delta p^{\mu}_{j^{\prime}}$
$\displaystyle=$
$\displaystyle\sum_{\mu\in{\mathcal{M}}}\sum_{{j^{\prime}}\in{\mathcal{P}}}\frac{\partial\hat{v}_{i}(K,{\bf
I})}{\partial p^{\mu}_{j^{\prime}}}P_{j{j^{\prime}}}\Delta S_{\mu j}.$

From Eq. (S7) and $\hat{{{\bf v}}}=P\hat{{{\bf f}}}$,

$\displaystyle\Delta\hat{v}_{i}(\\{\Delta S_{\alpha
j}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}})$ $\displaystyle=$
$\displaystyle-\sum_{\mu\in{\mathcal{M}}}\sum_{{j^{\prime}}\in{\mathcal{P}}}\Delta
S_{\mu
j}P_{j{j^{\prime}}}\hat{f}_{{j^{\prime}}}\frac{\partial\hat{v}_{i}}{\partial
I_{\mu}}.$ (S8) $\displaystyle=$
$\displaystyle-\hat{v}_{j}\sum_{\mu\in{\mathcal{M}}}\Delta S_{\mu
j}\frac{\partial\hat{v}_{i}}{\partial I_{\mu}}.$

Then, by defining $c^{\nu}_{\mu}(i):=\frac{\partial\hat{v}_{i}}{\partial
I_{\mu}}/\frac{\partial\hat{v}_{i}}{\partial I_{\nu}}$ and

$\displaystyle\Delta q_{i}^{\nu}(\\{\Delta S_{\mu
i}\\}_{\mu\in{\mathcal{M}}}):=\sum_{\mu\in{\mathcal{M}}}c^{\nu}_{\mu}(i)\Delta
S_{\mu i}$

and applying it to Eq. (S8) in the case of $j=i$, we obtain

$\displaystyle\Delta\hat{v}_{i}(\\{\Delta S_{\alpha
j}\\}_{\alpha\in{\mathcal{M}}\cup{\mathcal{C}}})=-\hat{v}_{i}\frac{\partial\hat{v}_{i}}{\partial
I_{\nu}}\underbrace{\sum_{\mu\in{\mathcal{M}}}c^{\nu}_{\mu}(i)\Delta S_{\mu
i}}_{=:\Delta q_{i}^{\nu}}.$

This is the derivation of Eqs. (4-5) in the main text.

Figure S2:  Numerical calculation with the genome-scale E.coli iJR904 model
reed2003expanded  for responses of the optimized reaction fluxes $\hat{{{\bf
v}}}$. The metabolic responses to effective changes in the input
stoichiometric coefficients $\\{S_{\mu i}\\}_{\mu\in{\mathcal{M}}}$,
$\Delta\hat{v}_{i}(\\{\Delta S_{\mu i}\\}_{\mu\in{\mathcal{M}}})/\Delta
q^{\mathrm{Glc}}_{i}$, are plotted against those to intake changes,
$\hat{v}_{i}\Delta\hat{v}_{i}(\Delta I_{\mathrm{Glc}})/\Delta
I_{\mathrm{Glc}}$. $I_{\mathrm{Glc}}\,{=}\,8.075\,\mathrm{[mmol/gDW/h]}$. Each
marker denotes a different reaction $i$.  Figure S3:  The metabolic responses
with several approximated estimations of $c^{\nu}_{\mu}$ are plotted. (Top)
With a constraint (i.e., the number of elements in ${\mathcal{C}}$,
$|{\mathcal{C}}|$, is 1). (Bottom) Without constraints (i.e.,
$|{\mathcal{C}}|$ is 0). From the left, $c^{\nu}_{\mu}$ is approximated as $1$
independent of $\mu$, $c^{\nu}_{\mu}$ is approximated as the ratios of the
carbon numbers of species $\mu$ and $\nu$, $c^{\nu}_{\mu}$ is approximated as
$\frac{\partial\Lambda}{\partial I_{\mu}}/\frac{\partial\Lambda}{\partial
I_{\nu}}$ (that equals to the ratio of the shadow price reznik2013flux  of
species $\mu$ and $\nu$, $y_{\mu}/y_{\nu}$), and Eq. (5) in the main text
(without approximation).  Figure S4:  The metabolic responses with matter
production as the objective reaction. $|{\mathcal{C}}|=0$. (A) Ethanol
production (acetaldehyde dehydrogenase; abbreviated as “ADHEr” in E. coli core
model). (B) ATP production (ATP maintenance requirement; abbreviated as “ATPM”
in E. coli core model).  Table 1: Description of the symbols in the main text

Symbol | Description
---|---
${\mathcal{M}},{\mathcal{C}},{\mathcal{R}},{\mathcal{P}}$ | Set of chemical species / constraints / reactions / pathways
$\mathcal{E},{\mathcal{O}}$ | Set of exchangeable species (${\mathcal{E}}\subset{\mathcal{M}}$) / objective components (${\mathcal{O}}\subset{\mathcal{M}}\cup{\mathcal{C}}$)
$S,K$ | Stoichiometry matrix for reactions / pathways ($S_{\alpha i},K_{\alpha{i^{\prime}}}\in\mathbb{R};\alpha\in{\mathcal{M}}\cup{\mathcal{C}},i\in{\mathcal{R}},{i^{\prime}}\in{\mathcal{P}}$)
$v_{i},f_{i^{\prime}}$ | Non-negative flux of reaction $i\,(\in{\mathcal{R}})$ / pathway $i^{\prime}\,(\in{\mathcal{P}})$
$q^{\nu}_{i},p^{\nu}_{i^{\prime}}$ | Metabolic price of reaction $i$ / pathway $i^{\prime}$ in terms of $\nu\,(\in{\mathcal{M}}\cup{\mathcal{C}}$)
$I_{\alpha}$ | Maximal intake of species $\alpha\,(\in{\mathcal{M}}$) or total capacity for constraint $\alpha\,(\in{\mathcal{C}}$)
$\Lambda$ | Flux of objective reaction $o\,(\in{\mathcal{R}})$