of probability landscapes when the estimated errors are sufficiently small, as
are the cases for the three examples given here. Overall, the goal of this
study is to provide a methodology for high precision solutions to the dCME
that can be applied to a large class of problems. Important unknown features
such as basins, attractors, and transitions for many biological networks can
be uncovered, analyzed, and their biological significance assessed. It is now
possible to analyze details of the topological, topographical, as well as
dynamic properties of the probability landscape for a large number of
biologically important stochastic networks that are previously not amenable to
computational investigations.

##  5 ACKNOWLEDGMENTS

This work was supported by NIH grant GM079804, NSF grant MCB1415589, and the
Chicago Biomedical Consortium with support from the Searle Funds at The
Chicago Community Trust. We thank Dr. Ao Ma for helpful comments, and Alan
Perez-Rathke for helpful discussions. YC acknowledges the support of the U.S.
Department of Energy through the LANL/LDRD Program for this work. YC thanks
Dr. Alan Perelson for helpful comments.

## References

  * [1] Ping Ao.  Laws in darwinian evolutionary theory.  Physics of life Reviews, 2(2):117–156, 2005.
  * [2] Jacob Stewart-Ornstein and Hana El-Samad.  Stochastic modeling of cellular networks.  Computational Methods in Cell Biology, 110:111, 2012.
  * [3] Hong Qian.  Cooperativity in cellular biochemical processes: noise-enhanced sensitivity, fluctuating enzyme, bistability with nonlinear feedback, and other mechanisms for sigmoidal responses.  Annual Review of Biophysics, 41:179–204, 2012.
  * [4] Adam Arkin, John Ross, and Harley H. McAdams.  Stochastic kinetic analysis of developmental pathway bifurcation in phage $\lambda$-infected Escherichia coli cells.  Genetics, 149(4):1633–1648, 1998.
  * [5] P.S. Swain, M.B. Elowitz, and E.D. Siggia.  Intrinsic and extrinsic contributions to stochasticity in gene expression.  Proceedings of the National Academy of Sciences of the United States of America, 99(20):12795, 2002.
  * [6] M.B. Elowitz, A.J. Levine, E.D. Siggia, and P.S. Swain.  Stochastic gene expression in a single cell.  Science, 297(5584):1183, 2002.
  * [7] Youfang Cao, Hsiao-Mei Lu, and Jie Liang.  Probability landscape of heritable and robust epigenetic state of lysogeny in phage lambda.  Proceedings of the National Academy of Sciences of the United States of America, 107(43):18445–18450, 2010.
  * [8] H.H. McAdams and A. Arkin.  It’s a noisy business! Genetic regulation at the nanomolar scale.  Trends in Genetics, 15(2):65–69, 1999.
  * [9] Darren J Wilkinson.  Stochastic modelling for quantitative description of heterogeneous biological systems.  Nature Reviews Genetics, 10(2):122–133, 2009.
  * [10] D. T. Gillespie.  Exact stochastic simulation of coupled chemical reactions.  Journal of Physical Chemistry, 81:2340–2361, 1977.
  * [11] Daniel T. Gillespie.  A rigorous derivation of the chemical master equation.  Physica A, 188:404–425, 1992.
  * [12] N.G. Van Kampen.  Stochastic processes in physics and chemistry, 3rd Edition.  Elsevier Science and Technology books, 2007.
  * [13] D.A. Beard and H. Qian.  Chemical biophysics: quantitative analysis of cellular systems.  Cambridge Univ Pr, 2008.
  * [14] Daniel T. Gillespie.  A diffusional bimolecular propensity function.  The Journal of Chemical Physics, 131(16):–, 2009.
  * [15] I.G. Darvey, B.W. Ninham, and P.J. Staff.  Stochastic models for second order chemical reaction kinetics. the equilibrium state.  The Journal of Chemical Physics, 45:2145–2155, 1966.
  * [16] D.A. McQuarrie.  Stochastic approach to chemical kinetics.  Journal of Applied Probability, 4:413–478, 1967.
  * [17] I.J. Laurenzi.  An analytical solution of the stochastic master equation for reversible bimolecular reaction kinetics.  The Journal of Chemical Physics, 113:3315–3322, 2000.
  * [18] Melissa Vellela and Hong Qian.  A quasistationary analysis of a stochastic chemical reaction: Keizer’s paradox.  Bulletin of Mathematical Biology, 69(5):1727–1746, 2007.
  * [19] N. G. Van Kampen.  A power series expansion of the master equation.  Canadian Journal of Physics, 39(4):551–567, 1961.
  * [20] Daniel. T. Gillespie.  The chemical langevin equation.  The Journal of Chemical Physics, 113:297–306, 2000.
  * [21] Daniel T. Gillespie.  The chemical Langevin and Fokker−Planck equations for the reversible isomerization reaction†.  The Journal of Physical Chemistry A, 106(20):5063–5071, 2002.
  * [22] Eric L. Haseltine and James B. Rawlings.  Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics.  The Journal of Chemical Physics, 117(15):6959–6969, 2002.
  * [23] C. W. Gardiner.  Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences.  Springer, New York, 2004.
  * [24] P Ao, D Galas, L Hood, L Yin, and XM Zhu.  Towards predictive stochastic dynamical modeling of cancer genesis and progression.  Interdisciplinary Sciences: Computational Life Sciences, 2(2):140–144, 2010.
  * [25] Jianghong Shi, Tianqi Chen, Ruoshi Yuan, Bo Yuan, and Ping Ao.  Relation of a new interpretation of stochastic differential equations to ito process.  Journal of Statistical Physics, 148(3):579–590, 2012.
  * [26] R. Grima, P. Thomas, and A. V. Straube.  How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations?  The Journal of Chemical Physics, 135(8):084103, Aug 2011.
  * [27] P. Thomas, H. Matuschek, and R. Grima.  How reliable is the linear noise approximation of gene regulatory networks?  BMC Genomics, 14 Suppl 4:S5, 2013.
  * [28] R.J. Allen, P.B. Warren, and P.R. Ten Wolde.  Sampling rare switching events in biochemical networks.  Physical Review Letters, 94(1):18104, 2005.
  * [29] H. Kuwahara and I. Mura.  An efficient and exact stochastic simulation method to analyze rare events in biochemical systems.  The Journal of Chemical Physics, 129:165101, 2008.
  * [30] B.J. Daigle, M.K. Roh, D.T. Gillespie, and L.R. Petzold.  Automated estimation of rare event probabilities in biochemical systems.  The Journal of Chemical Physics, 134:044110, 2011.
  * [31] Shuyun Jiao, Yanbo Wang, Bo Yuan, and Ping Ao.  Kinetics of muller’s ratchet from adaptive landscape viewpoint.  In Systems Biology (ISB), 2011 IEEE International Conference on, pages 27–32. IEEE, 2011.
  * [32] Youfang Cao and Jie Liang.  Adaptively biased sequential importance sampling for rare events in reaction networks with comparison to exact solutions from finite buffer dCME method.  The Journal of Chemical Physics, 139(2):025101, 2013.
  * [33] Gaowei Wang, Xiaomei Zhu, Jianren Gu, and Ping Ao.  Quantitative implementation of the endogenous molecular–cellular network hypothesis in hepatocellular carcinoma.  Interface focus, 4(3):20130064, 2014.
  * [34] Brian Munsky and Mustafa Khammash.  The finite state projection algorithm for the solution of the chemical master equation.  The Journal of Chemical Physics, 124(4):044104, 2006.
  * [35] Youfang Cao and Jie Liang.  Optimal enumeration of state space of finitely buffered stochastic molecular networks and exact computation of steady state landscape probability.  BMC Systems Biology, 2(1):30, 2008.
  * [36] Verena Wolf, Rushil Goel, Maria Mateescu, and Thomas Henzinger.  Solving the chemical master equation using sliding windows.  BMC Systems Biology, 4(1):42, 2010.
  * [37] Tobias Jahnke.  On reduced models for the chemical master equation.  Multiscale Modeling & Simulation, 9(4):1646–1676, 2011.
  * [38] Roger B Sidje.  Expokit: a software package for computing matrix exponentials.  ACM Transactions on Mathematical Software (TOMS), 24(1):130–156, 1998.
  * [39] Brian Munsky and Mustafa Khammash.  A multiple time interval finite state projection algorithm for the solution to the chemical master equation.  Journal of Computational Physics, 226(1):818 – 835, 2007.
  * [40] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, Clifford Stein, et al.  Introduction to algorithms, volume 2.  MIT press Cambridge, 2001.
  * [41] Fan RK Chung.  Spectral graph theory, volume 92.  American Mathematical Soc., 1997.
  * [42] Youfang Cao, Anna Terebus, and Jie Liang.  State space truncation with quantified errors for accurate solutions to discrete chemical master equation.  Bulletin of Mathematical Biology, page in press, 2016.
  * [43] H.M. Taylor and S. Karlin.  An Introduction to Stochastic Modeling, 3rd Ed. Academic Press, 1998.
  * [44] Bennett L Fox and Peter W Glynn.  Computing poisson probabilities.  Communications of the ACM, 31(4):440–445, 1988.
  * [45] Laurent Truffet.  Near complete decomposability: bounding the error by a stochastic comparison method.  Advances in Applied Probability, pages 830–855, 1997.
  * [46] A Irle.  Stochastic ordering for continuous-time processes.  Journal of Applied Probability, pages 361–375, 2003.
  * [47] Vladimir Kazeev, Mustafa Khammash, Michael Nip, and Christoph Schwab.  Direct solution of the chemical master equation using quantized tensor trains.  PLoS Computational Biology, 10(3):e1003359, 03 2014.
  * [48] R. Lehoucq, D. Sorensen, and C. Yang.  Arpack users’ guide: Solution of large scale eigenvalue problems with implicitly restarted Arnoldi methods.  SIAM, Philadelphia, 1998.
  * [49] W.J. Stewart.  Introduction to the numerical solution of Markov chains.  Princeton University Press NJ, 1994.
  * [50] Youcef Saad and Martin H Schultz.  GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986.
  * [51] Donald Gross and Douglas R. Miller.  The randomization technique as a modeling tool and solution procedure for transient Markov processes.  Operations Research, 32(2):343–361, 1984.
  * [52] Dan T Gillespie, Min Roh, and Linda R Petzold.  Refining the weighted stochastic simulation algorithm.  The Journal of Chemical Physics, 130(17):174103, 2009.
  * [53] Rory M Donovan, Andrew J Sedgewick, James R Faeder, and Daniel M Zuckerman.  Efficient stochastic simulation of chemical kinetics networks using a weighted ensemble of trajectories.  The Journal of Chemical Physics, 139(11):115105, 2013.
  * [54] Timothy S Gardner, Charles R Cantor, and James J Collins.  Construction of a genetic toggle switch in Escherichia coli.  Nature, 403(6767):339–342, 2000.
  * [55] Thomas B Kepler and Timothy C Elston.  Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations.  Biophysical Journal, 81(6):3116–3136, 2001.
  * [56] Keun-Young Kim and Jin Wang.  Potential energy landscape and robustness of a gene regulatory network: toggle switch.  PLoS Computational Biology, 3(3):e60, 2007.
  * [57] D. Schultz, J. N. Onuchic, and P. G. Wolynes.  Understanding stochastic simulations of the smallest genetic networks.  The Journal of Chemical Physics, 126(24):245102, 2007.
  * [58] Brian Munsky and Mustafa Khammash.  The finite state projection approach for the analysis of stochastic noise in gene networks.  Automatic Control, IEEE Transactions on, 53(Special Issue):201–214, 2008.
  * [59] Peter Deuflhard, Wilhelm Huisinga, T Jahnke, and Michael Wulkow.  Adaptive discrete Galerkin methods applied to the chemical master equation.  SIAM Journal on Scientific Computing, 30(6):2990–3011, 2008.
  * [60] Paul Sjöberg, Per Lötstedt, and Johan Elf.  Fokker–Planck approximation of the master equation in molecular biology.  Computing and Visualization in Science, 12(1):37–50, 2009.
  * [61] Moises Santillán.  On the use of the hill functions in mathematical models of gene regulatory networks.  Mathematical Modelling of Natural Phenomena, 3(02):85–97, 2008\.
  * [62] Haseong Kim and Erol Gelenbe.  Stochastic gene expression modeling with hill function for switch-like gene responses.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 9(4):973–979, 2012.
  * [63] M.K. Roh, B.J. Daigle, D.T. Gillespie, and L.R. Petzold.  State-dependent doubly weighted stochastic simulation algorithm for automatic characterization of stochastic biochemical rare events.  Journal of Chemical Physics, 135(23):234108, 2011.
  * [64] Mark Ptashne.  A Genetic Switch: Phage Lambda Revisited. Cold Spring Harbor Laboratory Press; 3 edition, 2004.
  * [65] Erik Aurell, Stanley Brown, Johan Johanson, and Kim Sneppen.  Stability puzzles in phage $\lambda$.  Physical Review E, 65(5):051914, 2002.
  * [66] Erik Aurell and Kim Sneppen.  Epigenetics as a first exit problem.  Physical Review Letters, 88(4):048101, 2002.
  * [67] X.-M. Zhu, L. Yin, L. Hood, and P. Ao.  Robustness, stability and efficiency of phage lambda genetic switch: dynamical structure analysis.  Journal of Bioinformatics and Computational Biology, 2:785–817, 2004.
  * [68] X.-M. Zhu, L. Yin, L. Hood, and P. Ao.  Calculating biological behaviors of epigenetic states in the phage $\lambda$ life cycle.  Functional & Integrative Genomics, 4(3):188–195, 2004.
  * [69] D. Hawley and W. McClure.  In vitro comparison of initiation properties of bacteriophage lambda wild-type PR and x3 mutant promoters.  Proceedings of the National Academy of Sciences of the United States of America, 77(11):6381–6385, 1980.
  * [70] D. Hawley and W. McClure.  Mechanism of activation of transcription initiation from the lambda PRM promoter.  Journal of Molecular Biology, 157(3):493–525, 1982.
  * [71] Madeline A. Shea and Gary K. Ackers.  The $OR$ control system of bacteriophage lambda a physical-chemical model for gene regulation.  Journal of Molecular Biology, 181(2):211–230, 1985.
  * [72] M. Li, W. McClure, and M. Susskind.  Changing the mechanism of transcriptional activation by phage lambda repressor.  Proceedings of the National Academy of Sciences of the United States of America, 94(8):3691–3696, 1997.
  * [73] Céline Kuttler and Joachim Niehren.  Gene Regulation in the Pi Calculus: Simulating Cooperativity at the Lambda Switch.  Transactions on Computational Systems Biology VII, 4230:24–55, 2006\.
  * [74] Gary L Johnson and Razvan Lapadat.  Mitogen-activated protein kinase pathways mediated by erk, jnk, and p38 protein kinases.  Science, 298(5600):1911–1912, 2002.
  * [75] Chi-Ying Huang and James E Ferrell.  Ultrasensitivity in the mitogen-activated protein kinase cascade.  Proceedings of the National Academy of Sciences of the United States of America, 93(19):10078–10083, 1996.
  * [76] Xiao Wang, Nan Hao, Henrik G Dohlman, and Timothy C Elston.  Bistability, stochasticity, and oscillations in the mitogen-activated protein kinase cascade.  Biophysical Journal, 90(6):1961–1978, 2006.
  * [77] Liang Qiao, Robert B Nachbar, Ioannis G Kevrekidis, and Stanislav Y Shvartsman.  Bistability and oscillations in the Huang-Ferrell model of MAPK signaling.  PLoS Computational Biology, 3(9):e184, 2007.
  * [78] N. I. Markevich, J. B. Hoek, and B. N. Kholodenko.  Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades.  The Journal of Cell Biology, 164(3):353–359, 2004.
  * [79] Paul Smolen, Douglas A. Baxter, and John H. Byrne.  Bistable MAP kinase activity: a plausible mechanism contributing to maintenance of late long-term potentiation.  American Journal of Physiology - Cell Physiology, 294(2):C503–C515, 2008.
  * [80] Marco Maggioni, Tanya Berger-Wolf, and Jie Liang.  Gpu-based steady-state solution of the chemical master equation.  In Parallel and Distributed Processing Symposium Workshops & PhD Forum (IPDPSW), 2013 IEEE 27th International, pages 579–588. IEEE, 2013\.
  * [81] John W. Little, Donald P. Shepley, and David W. Wert.  Robustness of a gene regulatory circuit.  The EMBO Journal, 18(15):4299–4307, 1999.
  * [82] Douglas Hanahan and Robert A Weinberg.  The hallmarks of cancer.  Cell, 100(1):57–70, 2000.
  * [83] Rosalind J Allen, Chantal Valeriani, and Pieter Rein ten Wolde.  Forward flux sampling for rare event simulations.  Journal of Physics: Condensed Matter, 21(46):463102, 2009.
  * [84] Russell Merris.  Laplacian matrices of graphs: a survey.  Linear Algebra and Its Applications, 197:143–176, 1994.

## Appendix

### Graph of Reaction Network, its Adjacency and Laplacian Matrices

${\mbox{\boldmath$G$}}_{R}$ can be represented by an $m\times m$ adjacency
matrix $C$, where:

${\mbox{\boldmath$C$}}^{m\times
m}=||{\mbox{\boldmath$C$}}_{i,\,j}||=\left\\{\begin{array}[]{l
l}1,&\quad\text{if $e_{ij}$ exists,}\\\ 0,&\quad\text{otherwise}\\\
\end{array}\right.$ (28)

The diagonal degree matrix $D$ of the graph ${\mbox{\boldmath$G$}}_{R}$ is:

${\mbox{\boldmath$D$}}^{m\times
m}=||{\mbox{\boldmath$D$}}_{i,\,j}||=\left\\{\begin{array}[]{l
l}\sum_{k=1}^{m}{{\mbox{\boldmath$C$}}_{i,k}},&\quad\text{if $i=j$},\\\
0,&\quad\text{if $i\neq j$},\\\ \end{array}\right.$ (29)

where each diagonal element ${\mbox{\boldmath$D$}}_{i,i}$ is the vertex degree
of the corresponding reaction $R_{i}$. The Laplacian matrix $L$ of the graph
${\mbox{\boldmath$G$}}_{R}$ can be then written as [84]:

${\mbox{\boldmath$L$}}={\mbox{\boldmath$D$}}-{\mbox{\boldmath$C$}}.$ (30)

Table 2: Detailed reactions and rate constants of genetic toggle switch. $R_{1}:GeneX\stackrel{{\scriptstyle k_{1}}}{{\rightarrow}}GeneX+X,\quad k_{1}=50\,s^{-1}$ | $R_{3}:X\stackrel{{\scriptstyle k_{3}}}{{\rightarrow}}\emptyset,\quad k_{3}=1\,s^{-1}$
---|---
$R_{2}:GeneY\stackrel{{\scriptstyle k_{2}}}{{\rightarrow}}GeneY+Y,\quad k_{2}=100\,s^{-1}$ | $R_{4}:Y\stackrel{{\scriptstyle k_{4}}}{{\rightarrow}}\emptyset,\quad k_{4}=1\,s^{-1}$
$R_{5}:2X+GeneY\stackrel{{\scriptstyle k_{5}}}{{\rightarrow}}BGeneY,\quad k_{5}=1\times 10^{-5}\,nM^{-2}\cdot s^{-1}$ | $R_{7}:BGeneY\stackrel{{\scriptstyle k_{7}}}{{\rightarrow}}2X+GeneY,\quad k_{7}=0.1\,s^{-1}$
$R_{6}:2Y+GeneX\stackrel{{\scriptstyle k_{6}}}{{\rightarrow}}BGeneX,\quad k_{6}=1\times 10^{-5}\,nM^{-2}\cdot s^{-1}$ | $R_{8}:BGeneX\stackrel{{\scriptstyle k_{8}}}{{\rightarrow}}2Y+GeneX,\quad k_{8}=0.1\,s^{-1}$
Table 3: Detailed reactions and rate constants of phage lambda epigenetic switch. $R_{1}:\emptyset+(OR3+OR2)\stackrel{{\scriptstyle s_{CI}^{0}}}{{\rightarrow}}CI+(OR3+OR2),s_{CI}^{0}=0.0069/s,$ | $R_{7}:2CI+OR1\stackrel{{\scriptstyle b_{CI}}}{{\rightarrow}}ROR1,b_{CI}=0.0021/nM^{2}\cdot s,$
---|---
$R_{2}:\emptyset+(OR3+COR2)\stackrel{{\scriptstyle s_{CI}^{0}}}{{\rightarrow}}CI+(OR3+COR2),s_{CI}^{0}=0.0069/s,$ | $R_{8}:2CI+OR2\stackrel{{\scriptstyle b_{CI}}}{{\rightarrow}}ROR2,b_{CI}=0.0021/nM^{2}\cdot s,$
$R_{3}:\emptyset+(OR3+ROR2)\stackrel{{\scriptstyle s_{CI}^{1}}}{{\rightarrow}}CI+(OR3+ROR2),s_{CI}^{1}=0.066/s,$ | $R_{9}:2CI+OR3\stackrel{{\scriptstyle b_{CI}}}{{\rightarrow}}ROR3,b_{CI}=0.0021/nM^{2}\cdot s,$
$R_{4}:\emptyset+(OR1+OR2)\stackrel{{\scriptstyle s_{Cro}}}{{\rightarrow}}Cro+(OR1+OR2),s_{Cro}=0.0929/s,$ | $R_{10}:2Cro+OR1\stackrel{{\scriptstyle b_{Cro}}}{{\rightarrow}}COR1,b_{Cro}=0.01289/nM^{2}\cdot s,$
$R_{5}:CI\stackrel{{\scriptstyle d_{CI}}}{{\rightarrow}}\emptyset,d_{CI}=0.0027/s,$ | $R_{11}:2Cro+OR2\stackrel{{\scriptstyle b_{Cro}}}{{\rightarrow}}COR2,b_{Cro}=0.01289/nM^{2}\cdot s,$
$R_{6}:Cro\stackrel{{\scriptstyle d_{Cro}}}{{\rightarrow}}\emptyset,d_{Cro}=0.0025/s,$ | $R_{12}:2Cro+OR3\stackrel{{\scriptstyle b_{Cro}}}{{\rightarrow}}COR3,b_{Cro}=0.01289/nM^{2}\cdot s,$
$R_{13}:ROR1+(OR2)\stackrel{{\scriptstyle
u_{ROR1}}}{{\rightarrow}}2CI+OR1+(OR2),u_{ROR1}=0.03998/s,$
$R_{14}:ROR1+(ROR2+OR3)\stackrel{{\scriptstyle
u_{ROR1}^{12}}}{{\rightarrow}}2CI+OR1+(ROR2+OR3),u_{ROR1}^{12}=0.0005/s,$
$R_{15}:ROR1+(ROR2+ROR3)\stackrel{{\scriptstyle
u_{ROR1}^{123}}}{{\rightarrow}}2CI+OR1+(ROR2+ROR3),u_{ROR1}^{123}=0.05531/s,$
$R_{16}:ROR1+(ROR2+COR3)\stackrel{{\scriptstyle
u_{ROR1}^{12}}}{{\rightarrow}}2CI+OR1+(ROR2+COR3),u_{ROR1}^{12}=0.0005/s,$
$R_{17}:ROR1+(COR2)\stackrel{{\scriptstyle
u_{ROR1}}}{{\rightarrow}}2CI+OR1+(COR2),u_{ROR1}=0.03998/s,$
$R_{18}:ROR2+(OR1+OR3)\stackrel{{\scriptstyle
u_{ROR2}}}{{\rightarrow}}2CI+OR2+(OR1+OR3),u_{ROR2}=1.026/s,$
$R_{19}:ROR2+(ROR1+OR3)\stackrel{{\scriptstyle
u_{ROR2}^{12}}}{{\rightarrow}}2CI+OR2+(ROR1+OR3),u_{ROR2}^{12}=0.01284/s,$
$R_{20}:ROR2+(OR1+ROR3)\stackrel{{\scriptstyle
u_{ROR2}^{23}}}{{\rightarrow}}2CI+OR2+(OR1+ROR3),u_{ROR2}^{23}=0.00928/s,$
$R_{21}:ROR2+(ROR1+ROR3)\stackrel{{\scriptstyle
u_{ROR2}^{123}}}{{\rightarrow}}2CI+OR2+(ROR1+ROR3),u_{ROR2}^{123}=0.01284/s,$
$R_{22}:ROR2+(COR1+OR3)\stackrel{{\scriptstyle
u_{ROR2}}}{{\rightarrow}}2CI+OR2+(COR1+OR3),u_{ROR2}=1.026/s,$
$R_{23}:ROR2+(OR1+COR3)\stackrel{{\scriptstyle
u_{ROR2}}}{{\rightarrow}}2CI+OR2+(OR1+COR3),u_{ROR2}=1.026/s,$
$R_{24}:ROR2+(COR1+COR3)\stackrel{{\scriptstyle
u_{ROR2}}}{{\rightarrow}}2CI+OR2+(COR1+COR3),u_{ROR2}=1.026/s,$
$R_{25}:ROR2+(ROR1+COR3)\stackrel{{\scriptstyle
u_{ROR2}^{12}}}{{\rightarrow}}2CI+OR2+(ROR1+COR3),u_{ROR2}^{12}=0.01284/s,$
$R_{26}:ROR2+(COR1+ROR3)\stackrel{{\scriptstyle
u_{ROR2}^{23}}}{{\rightarrow}}2CI+OR2+(COR1+ROR3),u_{ROR2}^{23}=0.00928/s,$
$R_{27}:ROR3+(OR2)\stackrel{{\scriptstyle
u_{ROR3}}}{{\rightarrow}}2CI+OR3+(OR2),u_{ROR3}=5.19753/s,$
$R_{28}:ROR3+(ROR2+OR1)\stackrel{{\scriptstyle
u_{ROR3}^{23}}}{{\rightarrow}}2CI+OR3+(ROR2+OR1),u_{ROR3}^{23}=0.04702/s,$
$R_{29}:ROR3+(ROR2+ROR1)\stackrel{{\scriptstyle
u_{ROR3}^{123}}}{{\rightarrow}}2CI+OR3+(ROR2+ROR1),u_{ROR3}^{123}=5.19753/s,$
$R_{30}:ROR3+(ROR2+COR1)\stackrel{{\scriptstyle
u_{ROR3}^{23}}}{{\rightarrow}}2CI+OR3+(ROR2+COR1),u_{ROR3}^{23}=0.04702/s,$
$R_{31}:ROR3+(COR2)\stackrel{{\scriptstyle
u_{ROR3}}}{{\rightarrow}}2CI+OR3+(COR2),u_{ROR3}=5.19753/s,$
$R_{32}:COR1+(OR2)\stackrel{{\scriptstyle
u_{COR1}}}{{\rightarrow}}2Cro+OR1+(OR2),u_{COR1}=0.08999/s,$
$R_{33}:COR1+(ROR2)\stackrel{{\scriptstyle
u_{COR1}}}{{\rightarrow}}2Cro+OR1+(ROR2),u_{COR1}=0.08999/s,$
$R_{34}:COR1+(COR2+OR3)\stackrel{{\scriptstyle
u_{COR1}^{12}}}{{\rightarrow}}2Cro+OR1+(COR2+OR3),u_{COR1}^{12}=0.01776/s,$
$R_{35}:COR1+(COR2+ROR3)\stackrel{{\scriptstyle
u_{COR1}^{12}}}{{\rightarrow}}2Cro+OR1+(COR2+ROR3),u_{COR1}^{12}=0.01776/s,$
$R_{36}:COR1+(COR2+COR3)\stackrel{{\scriptstyle
u_{COR1}^{123}}}{{\rightarrow}}2Cro+OR1+(COR2+COR3),u_{COR1}^{123}=0.05531/s,$
$R_{37}:COR2+(OR1+OR3)\stackrel{{\scriptstyle
u_{COR2}}}{{\rightarrow}}2Cro+OR2+(OR1+OR3),u_{COR2}=0.6306/s,$
$R_{38}:COR2+(ROR1+OR3)\stackrel{{\scriptstyle
u_{COR2}}}{{\rightarrow}}2Cro+OR2+(ROR1+OR3),u_{COR2}=0.6306/s,$
$R_{39}:COR2+(OR1+ROR3)\stackrel{{\scriptstyle
u_{COR2}}}{{\rightarrow}}2Cro+OR2+(OR1+ROR3),u_{COR2}=0.6306/s,$
$R_{40}:COR2+(ROR1+ROR3)\stackrel{{\scriptstyle
u_{COR2}}}{{\rightarrow}}2Cro+OR2+(ROR1+ROR3),u_{COR2}=0.6306/s,$
$R_{41}:COR2+(COR1+OR3)\stackrel{{\scriptstyle
u_{COR2}^{12}}}{{\rightarrow}}2Cro+OR2+(COR1+OR3),u_{COR2}^{12}=0.12448/s,$
$R_{42}:COR2+(OR1+COR3)\stackrel{{\scriptstyle
u_{COR2}^{23}}}{{\rightarrow}}2Cro+OR2+(OR1+COR3),u_{COR2}^{23}=0.23822/s,$
$R_{43}:COR2+(COR1+COR3)\stackrel{{\scriptstyle
u_{COR2}^{123}}}{{\rightarrow}}2Cro+OR2+(COR1+COR3),u_{COR2}^{123}=0.14641/s,$
$R_{44}:COR2+(ROR1+COR3)\stackrel{{\scriptstyle
u_{COR2}^{23}}}{{\rightarrow}}2Cro+OR2+(ROR1+COR3),u_{COR2}^{23}=0.23822/s,$
$R_{45}:COR2+(COR1+ROR3)\stackrel{{\scriptstyle
u_{COR2}^{12}}}{{\rightarrow}}2Cro+OR2+(COR1+ROR3),u_{COR2}^{12}=0.12448/s,$
$R_{46}:COR3+(OR2)\stackrel{{\scriptstyle
u_{COR3}}}{{\rightarrow}}2Cro+OR3+(OR2),u_{COR3}=0.00928/s,$
$R_{47}:COR3+(ROR2)\stackrel{{\scriptstyle
u_{COR3}}}{{\rightarrow}}2Cro+OR3+(ROR2),u_{COR3}=0.00928/s,$
$R_{48}:COR3+(COR2+OR1)\stackrel{{\scriptstyle
u_{COR3}^{23}}}{{\rightarrow}}2Cro+OR3+(COR2+OR1),u_{COR3}^{23}=0.00351/s,$
$R_{49}:COR3+(COR2+ROR1)\stackrel{{\scriptstyle
u_{COR3}^{23}}}{{\rightarrow}}2Cro+OR3+(COR2+ROR1),u_{COR3}^{23}=0.00351/s,$
$R_{50}:COR3+(COR2+COR1)\stackrel{{\scriptstyle
u_{COR3}^{123}}}{{\rightarrow}}2Cro+OR3+(COR2+COR1),u_{COR3}^{123}=0.01092/s$
Table 4: Molecular species in the network of bistable MAPK signaling cascade. Molecular species | Descriptions
---|---
MEK | ERK kinase
MKP3 | ERK phosphatase
K | ERK, extracellular signal-regulated kinase
KpY | Single phosphorylated ERK on Y residue
KpT | Single phosphorylated ERK on T residue
Kpp | Dual phosphorylated ERK on both Y and T residue
K_MEK_Y | K bound by MEK at residue Y
K_MEK_T | K bound by MEK at residue T
KpY_MEK | KpY bound by MEK
KpT_MEK | KpT bound by MEK
Kpp_MKP3 | Kpp associated with MKP3
KpY_MKP3 | KpY associated with MKP3
KpT_MKP3_Y | KpT associated with MKP3 at residue Y
KpT_MKP3_T | KpT associated with MKP3 at residue T
K_MKP3_T | K associated with MKP3 at residue T
K_MKP3_Y | K associated with MKP3 at residue Y
Figure 13: A simplified conceptual model of the MAPK network. The MEK and ERK (K) form a positive feedback loop.  Table 5: Detailed reactions and rate constants in MAPK signaling network. $R_{1}:\emptyset\overset{\,s_{2}\,}{\underset{\,d_{2}\,}{\rightleftharpoons}}\text{MEK},s_{2}=0.001/s,d_{2}=0.15/s,$ | $R_{2}:\emptyset+\text{(Kpp)}\stackrel{{\scriptstyle s_{3}}}{{\rightarrow}}\text{MEK}+\text{(Kpp)},s_{3}=0.005/s,$
---|---
$R_{5}:\text{K}+\text{MEK}\overset{\,k_{1}\,}{\underset{\,k_{-1}\,}{\rightleftharpoons}}\text{K\\_MEK\\_Y},k_{1}=0.375/nM\cdot s,k_{-1}=1.0/s,$ | $R_{3}:\emptyset\overset{\,s_{1}\,}{\underset{\,d_{1}\,}{\rightleftharpoons}}\text{K},s_{1}=0.00024/s,d_{1}=0.0001/s,$
$R_{7}:\text{KpY}+\text{MEK}\overset{\,k_{3}\,}{\underset{\,k_{-3}\,}{\rightleftharpoons}}\text{KpY\\_MEK},k_{3}=0.375/nM\cdot s,k_{-3}=1.0/s,$ | $R_{4}:\text{KpY}\stackrel{{\scriptstyle d_{1}}}{{\rightarrow}}\emptyset,\text{KpT}\stackrel{{\scriptstyle d_{1}}}{{\rightarrow}}\emptyset,\text{Kpp}\stackrel{{\scriptstyle d_{1}}}{{\rightarrow}}\emptyset,d_{1}=0.0001/s,$
$R_{9}:\text{K}+\text{MEK}\overset{\,k_{5}\,}{\underset{\,k_{-5}\,}{\rightleftharpoons}}\text{K\\_MEK\\_T},k_{5}=0.375/nM\cdot s,k_{-5}=1.0/s,$ | $R_{6}:\text{K\\_MEK\\_Y}\stackrel{{\scriptstyle k_{2}}}{{\rightarrow}}\text{KpY}+\text{MEK},k_{2}=0.06/s,$
$R_{11}:\text{KpT}+\text{MEK}\overset{\,k_{7}\,}{\underset{\,k_{-7}\,}{\rightleftharpoons}}\text{KpT\\_MEK},k_{7}=0.375/nM\cdot s,k_{-7}=1.0/s,$ | $R_{8}:\text{KpY\\_MEK}\stackrel{{\scriptstyle k_{4}}}{{\rightarrow}}\text{Kpp}+\text{MEK},k_{4}=4.5/s,$
$R_{13}:\text{Kpp}+\text{MKP3}\overset{\,h_{1}\,}{\underset{\,h_{-1}\,}{\rightleftharpoons}}\text{Kpp\\_MKP3},h_{1}=0.015/nM\cdot s,h_{-1}=1.0/s,$ | $R_{10}:\text{K\\_MEK\\_T}\stackrel{{\scriptstyle k_{6}}}{{\rightarrow}}\text{KpT}+\text{MEK},k_{6}=0.06/s,$
$R_{15}:\text{KpT\\_MKP3\\_Y}\overset{\,h_{3}\,}{\underset{\,h_{-3}\,}{\rightleftharpoons}}\text{KpT}+\text{MKP3},h_{3}=0.31/s,h_{-3}=0.01/nM\cdot s,$ | $R_{12}:\text{KpT\\_MEK}\stackrel{{\scriptstyle k_{8}}}{{\rightarrow}}\text{Kpp}+\text{MEK},k_{8}=4.5/s,$
$R_{16}:\text{KpT}+\text{MKP3}\overset{\,h_{4}\,}{\underset{\,h_{-4}\,}{\rightleftharpoons}}\text{KpT\\_MKP3\\_T},h_{4}=0.01/nM\cdot s,h_{-4}=1.0/s,$ | $R_{14}:\text{Kpp\\_MKP3}\stackrel{{\scriptstyle h_{2}}}{{\rightarrow}}\text{KpT\\_MKP3\\_Y},h_{2}=0.032/s,$
$R_{18}:\text{K\\_MKP3\\_T}\overset{\,h_{6}\,}{\underset{\,h_{-6}\,}{\rightleftharpoons}}\text{K}+\text{MKP3},h_{6}=0.086/s,h_{-6}=0.0011/nM\cdot s,$ | $R_{17}:\text{KpT\\_MKP3\\_T}\stackrel{{\scriptstyle h_{5}}}{{\rightarrow}}\text{K\\_MKP3\\_T},h_{5}=0.5/s,$
$R_{19}:\text{KpY}+\text{MKP3}\overset{\,h_{7}\,}{\underset{\,h_{-7}\,}{\rightleftharpoons}}\text{KpY\\_MKP3},h_{7}=0.01/nM\cdot s,h_{-7}=1.0/s,$ | $R_{20}:\text{KpY\\_MKP3}\stackrel{{\scriptstyle h_{8}}}{{\rightarrow}}\text{K\\_MKP3\\_Y},h_{8}=0.47/s,$
$R_{21}:\text{K\\_MKP3\\_Y}\overset{\,h_{9}\,}{\underset{\,h_{-9}\,}{\rightleftharpoons}}\text{K}+\text{MKP3},h_{9}=0.14/s,h_{-9}=0.0018/nM\cdot s.$ |


# PartonDensity.jl: a novel parton density determination code

Francesca Capel<EMAIL_ADDRESS>Max-Planck-Institut für Physik, München,
Germany  Ritu Aggarwal<EMAIL_ADDRESS>USAR, Guru Gobind Singh
Indraprastha University, East Delhi-110032  Michiel Botje<EMAIL_ADDRESS>Nikhef, Amsterdam, The Netherlands  Allen Caldwell<EMAIL_ADDRESS>Max-
Planck-Institut für Physik, München, Germany  Oliver Schulz
<EMAIL_ADDRESS>Max-Planck-Institut für Physik, München, Germany  Andrii
Verbytskyi<EMAIL_ADDRESS>Max-Planck-Institut für Physik,
München, Germany

###### Abstract

We introduce our novel Bayesian parton density determination code,
PartonDensity.jl. The motivation for this new code, the framework and its
validation are described. As we show, PartonDensity.jl provides both a
flexible environment for the determination of parton densities and a wealth of
information concerning the knowledge update provided by the analyzed data set.

##  I Introduction

This paper describes a novel approach to the determination of parton density
functions (PDFs) of hadrons and presents the PartonDensity.jl programming code
that implements this approach.

Our analysis method significantly differs from those pursued by other groups
Cridge  _et al._ (2022); Bailey  _et al._ (2021); Hou  _et al._ (2021); Dulat
_et al._ (2016); Alekhin  _et al._ (2018, 2014); Ball  _et al._ (2022, 2017);
Abt  _et al._ (2022); Abramowicz  _et al._ (2015), allowing us to tackle data
sets that have not been included in PDF determinations so far. The features
setting our code apart from existing ones are described in the following.

First of all, we use Bayesian techniques to determine the parton density
functions.111 For other efforts towards PDF determination from a Bayesian
perspective, see e.g. Cocuzza  _et al._ (2021); Hunt-Smith  _et al._ (2022).
This formulation allows us to include in a coherent and transparent way known
constraints and results from previous theoretical and experimental analyses in
prior probability distributions. The result of a Bayesian analysis is a
multivariate posterior probability distribution of all model parameters,
including nuisance parameters used to describe systematic uncertainties in the
data. Correlations between any subset of parameters can be studied, providing
vastly more information than results obtained otherwise. Systematic error
propagation is simply achieved by integrating the posterior over the nuisance
parameters. In addition, the information content of the data is easily judged
by comparing the posterior to the input priors.

Secondly, we have implemented a forward modeling approach, where event numbers
in kinematic bins are predicted and compared to observed event counts, which
are always bin-by-bin uncorrelated with known Poisson distribution. This is in
contrast to the analyses of unfolded cross sections, where a Gaussian
statistical hypothesis is implied and where the data are always correlated in
a way that is often not known. Another advantage of forward modeling is that
it can correctly handle the low statistics case with empty or sparsely
populated bins.

The current version of our code allows for the analysis of inclusive
$\mathrm{e}^{\pm}\mathrm{p}$ scattering data, made available in the form of
binned event counts. An extension to the analysis of differential cross
sections extracted at the QED-Born level is in development and will be
reported separately.

The PartonDensity.jl software package222 Available at
https://github.com/cescalara/PartonDensity.jl.  is implemented in the modern
julia Bezanson  _et al._ (2017) language and uses several julia packages for
the analysis, the most important of which is the Bayesian Analysis Toolkit
(BAT.jl) Schulz  _et al._ (2021). In the development of the code, a large
amount of effort was dedicated to speeding up the necessary calculations, and
we briefly report on these in this work.

This paper is organized as follows. In Section II, we begin with a general
introduction to the Bayesian analysis approach. In Section III, we describe
how the $\mathrm{e}^{\pm}\mathrm{p}$ binned event counts are computed,
starting from a set of parameterized input PDFs. Some technical details about
the implementation of our approach are given in Section IV. In Section V, we
discuss PDF parameterizations and prior parameter constraints. Finally,
validation tests on pseudo-data and goodness-of-fit tests are described in
Sections VI and VII, respectively. For this, we also refer the reader to
Aggarwal  _et al._ (2023) where we have published an analysis of the ZEUS data
Abramowicz  _et al._ (2014); Abt  _et al._ (2020) which very closely resemble
the pseudo-data used in the validation described here.

##  II Bayesian Analysis Approach

###  II.1 Introduction

We use a Bayesian approach to obtain the joint posterior distribution
$p(\theta|D,M)$ of a set of model parameters $\theta$, conditional to the data
set $D$ being analyzed. It is also conditional to the particular model choices
$M$ such as, among others, the modeling of systematic uncertainties in the
data and the specific form of the PDFs chosen.

Because it is not possible to obtain an analytic expression for the posterior,
Monte Carlo techniques are used to create parameter samples that are
distributed according to $p(\theta|D,M)$. Nuisance parameters that are needed
in the modeling of the data but are of no interest are removed by integrating
the posterior over these. In this way, systematic uncertainties are
automatically propagated in a consistent fashion by integrating over the
systematic parameters in the model.

The resulting samples then enable us to study single- and multi-dimensional
marginal distributions of the parameters of interest and extract all kinds of
statistical estimators like mean values, credibility intervals, and so forth.

The posterior parameter distribution is obtained from Bayes’ theorem,

$p(\theta|D,M)=\frac{p(D|\theta,M)\ p(\theta|M)}{p(D|M)}\;.$ (1)

Here $p(\theta|M)$ is the prior distribution of the parameters $\theta$ and
$p(D|\theta,M)$ is the probability to observe the data $D$, given a particular
set of parameter values $\theta$. In Eq. 1 it is evaluated as a function of
$\theta$ for fixed data $D$ and is called the likelihood function (it is not a
probability density), denoted as

$\mathcal{L}_{D}(\theta)=p(D|\theta),\quad\text{$D$ fixed}.$

Here and in the following, we will always imply the model choice $M$ but for
clarity, it will be dropped in the notation.

An important feature of Bayes’ theorem is the fact that the support of the
posterior can never exceed that of the prior. This makes it straightforward to
impose hard constraints on allowed parameter values, such as the requirement
that a parameter is positive definite.

The normalization of the posterior is given by the so-called evidence

$p(D)=\int\mathcal{L}_{D}(\theta)\ p(\theta)\ \,\mathrm{d}\theta\;.$ (2)

The evidence is a scalar value that is often challenging to compute as it
results from a multi-dimensional integral. However, as it is independent of
the parameter vector $\theta$, its value is not needed to study the relative
credibility of parameter values. Sampling algorithms like Markov Chain Monte
Carlo (MCMC) make it possible to draw samples distributed according to the
posterior distribution based on the non-normalized posterior density

$\tilde{p}(\theta|D)=\mathcal{L}_{D}(\theta)\ p(\theta)$ (3)

This yields posterior samples that can be used to produce the desired marginal
parameter distributions, which can always be normalized afterward.

###  II.2 Implementation of MCMC

In our experience, prior parameter distributions with hard constraints, like
positive definite parameters that may assume values close to zero can
negatively impact MCMC convergence and efficiency. Several of the parameters
of our model fall into this category. Furthermore, the momentum sum rule (see
below) introduces a strict correlation between the momentum parameters that
restricts the support of the posterior density to a subspace of the full
parameter space. Densities with such a support typically cause sampling
algorithms to fail completely, as the posterior density is zero almost
everywhere in the parameter space. The methods to circumvent this obstacle are
described below.

We solve the problem of sampling from a complicated distribution by performing
a change of variables with the aim of creating a parameter space that is
straightforward to sample from. Thus, instead of sampling the density
$\tilde{p}(\theta|D)$ directly, we introduce a suitable transformation
$\theta=f_{\mathrm{}}(x)$ and sample the density

$\tilde{p}_{\mathrm{}}(x|D)=\mathcal{L}_{D}[f_{\mathrm{}}(x)]\;\mathcal{N}(x)$
(4)

where $\mathcal{N}(x)$ is a multivariate Normal distribution of the
appropriate dimension. The posterior $\tilde{p}_{\mathrm{}}(x|D)$ is unbounded
and has full support over it’s parameter space, and so is much easier to
sample. Samples of the original parameters $\theta$ are then obtained by
simply applying the parameter transformation $f_{\mathrm{}}$ to generated
samples of the parameters $x$ in the alternate space.

The challenge lies in finding suitable transformations $f_{\mathrm{}}$. For
univariate components of the prior distribution, their quantile functions and
the cumulative distribution function of the normal distribution provide the
necessary building blocks Schulz  _et al._ (2021). For the Dirichlet
distribution that we use to satisfy the momentum sum rule in the prior (see
Section III.2), a suitable transformation is given in Betancourt (2012).

###  II.3 Marginalization and uncertainty propagation

In many model analyses one is interested not in the full posterior
distribution but in the marginal probability distribution of only one or a few
parameters. For example, the probability distribution of a single parameter
$\theta_{i}$ is obtained from

$p(\theta_{i}|D)=\int p(\theta|D)\,\mathrm{d}{\theta}_{j\neq i}\,.$

In practice, this calculation is performed quite simply by histogramming the
$\theta_{i}$ values from all posterior samples, ignoring the other parameters.
From such a histogram, several estimators can be computed to report the
results. Commonly used quantities are:

#### Mode of $\theta_{i}$.

The value ${\theta}_{i}^{*}$ where the marginalized posterior probability
density has a maximum. Modes are usually computed for fully marginalized
posteriors $p(\theta_{i}|D)$ (‘marginal mode’), for the entire posterior
(‘global mode’), or for any number of parameters, with the rest marginalized.
Note that the parameter values that maximize the full posterior distribution
usually do not coincide with those that maximize marginalized distributions.

#### Mean of $\theta_{i}$.

This is the expectation value

$<\theta_{i}>=\int_{\theta_{\mathrm{min}}}^{\theta_{\mathrm{max}}}p(\theta_{i}|D)\theta_{i}\,\mathrm{d}\theta_{i}\,,$

with the parameter bounds denoted by
$[\theta_{\mathrm{min}},\theta_{\mathrm{max}}]$.

#### Median of $\theta_{i}$.

The value $\hat{\theta}_{i}$ that splits the probability content of
$p(\theta_{i}|D)$ in two:

$\int_{\theta_{\mathrm{min}}}^{\hat{\theta}_{i}}p(\theta_{i}|D)\,\mathrm{d}\theta_{i}=\int_{\hat{\theta}_{i}}^{\theta_{\mathrm{max}}}p(\theta_{i}|D)\,\mathrm{d}\theta_{i}=0.5\,.$

#### Central Interval.

The $(1-2\alpha)$ central interval $[\theta_{\mathrm{-}},\theta_{\mathrm{+}}]$
is defined such that a fraction $\alpha$ of the probability is contained on
either side of the interval:

$\int_{\theta_{\mathrm{min}}}^{\theta_{\mathrm{-}}}p(\theta_{i}|D)\,\mathrm{d}\theta_{i}=\int_{\theta_{\mathrm{+}}}^{\theta_{\mathrm{max}}}p(\theta_{i}|D)\,\mathrm{d}\theta_{i}=\alpha.$

#### Smallest Interval.

The $\alpha$ smallest interval(s) is defined such that a fraction $\alpha$ of
the probability is contained in a set of intervals where the set size is
minimized. This is realized as a Lebesgue integral:

$\int_{p(\theta_{i}|D)\geq p_{\rm
min}}p(\theta_{i}|D)\,\mathrm{d}\theta_{i}=\alpha$

where $p_{\rm min}$ is to be determined. This procedure can result in several
intervals in $\theta_{i}$ in the case of multimodal distributions.

#### Uncertainty Propagation.

Having full access to the posterior allows for the evaluation of the
probability distribution of any function of the parameters. In contrast to
standard techniques used for error propagation, there is no need here for any
assumptions like distributions being Gaussian-shaped. As an example, consider
that we have a function $f(x|\theta)$ of interest, for instance, a parton
distribution that depends on a subset of the $\theta$ parameters. The
distribution of $f(x)$ is then given by

$\ p(f(x)|D)=\int p(f(x|\theta))p(\theta|D)\,\mathrm{d}\theta$

which can be evaluated in a straightforward way from the posterior samples. We
give examples of such uses below.

###  II.4 BAT.jl - the Bayesian Analysis Toolkit

For the Bayesian inference process, we use the software package BAT.jl Schulz
_et al._ (2021), which is a high-performance toolkit for Bayesian analysis
tasks, coded in the julia programming language Bezanson  _et al._ (2017). The
package provides multiple algorithms for posterior sampling, integration and
mode-finding, as well as automatic plotting and reporting functionality.

The BAT.jl package also has the ability to automatically generate suitable
space transformations between user-defined probability distributions and
standard multivariate normal or uniform distributions. This enables us to
automatically perform prior-space transformations as described in Section II.2
and sample the posterior via the Metropolis-Hastings algorithm in an
unconstrained space where the prior has become a normal distribution in all
dimensions. The samples are then automatically transformed back into the
original space.

A full description of the algorithms and tuning of the BAT.jl toolkit can be
found in Schulz  _et al._ (2021).

##  III Analysis Structure

As mentioned in the introduction, we focus in this initial release of our PDF
determination code on the analysis of electron(positron)-proton deep inelastic
scattering data. The data are generally reported in bins of the scaling
variables $x$ and $Q^{2}$, which are computed from the reconstructed four
momenta of the incoming proton and the incoming and scattered leptons (for
more details see, e.g., Abt  _et al._ (2020)). In this section, we describe
how the event numbers observed in these data are predicted, starting from a
set of proton PDFs parameterized at some input scale $Q^{2}_{0}$.

###  III.1 Computation of the $\mathrm{e}^{\pm}\mathrm{p}$ cross sections

To compute the $\mathrm{e}^{\pm}\mathrm{p}$ deep inelastic cross sections, we
start from a set of quark, antiquark, and gluon distributions $xf_{i}(x)$,
parameterized at a fixed value of $Q^{2}_{0}$. The aim of the analysis is to
determine from the data the posterior joint distribution of these parameters.

In the above, $f_{i}(x)$ is the number density of partons of type $i$ in the
proton, and $xf_{i}(x)$ is the momentum density of these partons.

The next step is to evolve these distributions in perturbative QCD Gribov and
Lipatov (1972a, b); Lipatov (1975); Dokshitzer (1977); Altarelli and Parisi
(1977) from the input scale to larger values of $Q^{2}$. We use the qcdnum
program Botje (2011) to evolve the PDFs in the $\overline{\mathrm{M}S}$-scheme
at next-to-next-to-leading order (NNLO) Tarasov  _et al._ (1980); Curci  _et
al._ (1980); Furmanski and Petronzio (1980, 1982); Larin and Vermaseren
(1993); Chetyrkin  _et al._ (1997); Buza  _et al._ (1998); Moch  _et al._
(2004); Vogt  _et al._ (2004); Ball  _et al._ (2016).

In the analysis, we impose the momentum sum rule and the valence quark
counting rules. The momentum sum rule states that the fractional momenta of
all partons in the proton add up to unity:

$\sum_{i}\int_{0}^{1}xf_{i}(x)\,\mathrm{d}x=\sum_{i}\Delta_{i}=1\;.$ (5)

We introduce here the notation $\Delta_{i}$ for the total momentum fraction
carried by the parton species $i$.

The quark counting rules fix the net number of quarks in the proton so that
its quantum numbers are conserved:

$\int_{0}^{1}[q_{i}(x)-\bar{q}_{i}(x)]\;\,\mathrm{d}x=\left\\{\begin{array}[]{l}2\
\ \text{for}\ \ i=u\;,\\\ 1\ \ \text{for}\ \ i=d\;,\\\ 0\ \ \text{for}\ \
i=s,c,b,t\;.\end{array}\right.$ (6)

Here and in the following we use the notation $q$, $\bar{q}$ and $g$ to denote
quark, antiquark, and gluon densities.

It is important to point out that the QCD evolution equations guarantee that
sum rules that are imposed at the starting scale $Q^{2}_{0}$ will be satisfied
at all scales.

The neutral current DIS cross section for $\mathrm{e}^{\pm}\mathrm{p}$
scattering is given in terms of generalized structure functions of the proton
as ($y$ is the inelasticity variable, see Abt  _et al._ (2020))

$\frac{{\,\mathrm{d}}^{2}\sigma^{\mathrm{e}^{\pm}\mathrm{p}}_{\mathrm{NC}}}{\,\mathrm{d}x\,\mathrm{d}Q^{2}}=\frac{2\pi\alpha}{xQ^{4}}(Y_{+}\tilde{F}_{2}\mp
Y_{-}x\tilde{F}_{3}-y^{2}\tilde{F}_{\mathrm{L}})$ (7)

where $Y_{\pm}=1\pm(1-y)^{2}$, and $\alpha$ is the fine-structure constant.

The generalized structure functions are related to the vector and axial-vector
coupling constants $v_{\mathrm{e}}$ and $a_{\mathrm{e}}$ by

$\begin{split}\tilde{F}_{i}&=F_{i}^{\gamma}-k_{\mathrm{Z}}v_{\mathrm{e}}\,F_{i}^{\gamma\mathrm{Z}}+k_{\mathrm{Z}}^{2}(v_{\mathrm{e}}^{2}+a_{\mathrm{e}}^{2})\,F_{i}^{\mathrm{Z}},\quad
i=2,\mathrm{L}\\\
x\tilde{F}_{3}&=k_{\mathrm{Z}}a_{\mathrm{e}}\,xF_{3}^{\gamma\mathrm{Z}}+k_{\mathrm{Z}}^{2}2v_{\mathrm{e}}a_{\mathrm{e}}\,xF_{3}^{\mathrm{Z}}\rule{0.0pt}{10.00002pt}\end{split}$
(8)

where

$k_{\mathrm{Z}}(Q^{2})=\frac{Q^{2}}{(Q^{2}+m_{\mathrm{Z}}^{2})\,4\sin^{2}{\theta_{\mathrm{w}}}\cos^{2}{\theta_{\mathrm{w}}}}.$

In our analysis, we use $\sin^{2}\theta_{\mathrm{w}}=0.231$ and
$m_{\mathrm{Z}}=91.1876$ GeV for the electroweak mixing angle and the Z-boson
mass, respectively Workman  _et al._ (2022).

In LO QCD, the structure functions are linear combinations of parton densities
of different flavor (note that $F_{\mathrm{L}}=0$ at LO).

$\begin{split}\\{F_{2}^{\gamma},F_{2}^{\gamma\mathrm{Z}},F_{2}^{Z}\\}&=x\sum_{i}\\{e_{i}^{2},2e_{i}g_{\mathrm{V}},g_{\mathrm{A}}^{2}\\}(q_{i}+\overline{q}_{i})\\\
x\\{F_{3}^{\gamma},F_{3}^{\gamma\mathrm{Z}},F_{3}^{Z}\\}&=x\sum_{i}\\{0,2e_{i}g_{A},2g_{\mathrm{V}}g_{\mathrm{A}}\\}(q_{i}-\overline{q}_{i})\end{split}$
(9)

where $e_{i}$ is the charge of the quark species $i$ and $g_{\mathrm{V}}$ and
$g_{\mathrm{A}}$ are the weak vector and axial-vector couplings of the quark
to the Z-boson.

We use the qcdnum package to compute the structure functions at NNLO which
involves convolutions of parton densities with various coefficient functions
Sanchez Guillen  _et al._ (1991); van Neerven and Zijlstra (1991); Zijlstra
and van Neerven (1991, 1992); van Neerven and Vogt (2000a, b). Note in this
respect that the coefficients in Eq. 9 are the same for all up-type (down-
type) quarks with $e_{i}=\frac{2}{3}\ (\frac{1}{3})$. Thus we can compute NNLO
structure functions separately for the sum of up-type and down-type quarks and
multiply these by the coefficients afterwards. This allows for an efficient
calculation of the cross section based on the NNLO computation of only 6
structure functions.

###  III.2 Forward model

In the forward modeling approach, we compute the expected number of events
$\nu$ in a bin $(\Delta x,\Delta Q^{2})$ as an integral of the differential
cross section over the full kinematic phase space $\mathcal{D}$. Introducing
the short-hand notation $\bm{u}=\\{x,Q^{2}\\}$ this integral can be written as

$\nu=\mathcal{L}\int_{\Delta\bm{u}^{\prime}}\left[\int_{\mathcal{D}}A(\bm{u}^{\prime}|\bm{u})\frac{\,\mathrm{d}^{2}\sigma(\bm{u})}{\,\mathrm{d}\bm{u}}\,\mathrm{d}\bm{u}\right]\,\mathrm{d}\bm{u^{\prime}}\,,$
(10)

where the primed (unprimed) variable refers to the reconstructed (true)
kinematics. Here $\mathcal{L}$ is the integrated luminosity of the data set
being analyzed, and $A$ is a transformation kernel that takes into account the
detector effects and radiative corrections to the QED Born-level differential
cross sections, computed as described in Section III.1. Note that the forward
model approach requires the experimental groups to provide the transformation
kernel as part of their data release.

We assume here that this information is cast in the form of a matrix $A$ that
provides a mapping from the Born-level cross section integrated over bins in
the kinematic domain to events accumulated in the experimental bins. Denoting
by $\nu_{i}$ the counts in an experimental bin $i$ and by $\lambda_{j}$ the
integrated cross section in a kinematic bin $j$, Eq. 10 is written as

$\nu_{i}=\mathcal{L}\sum_{j}A_{ij}\lambda_{j}$ (11)

Experimental systematic uncertainties are encoded in a set of deviation
matrices. The total systematic deviation is then given by a linear combination
of these deviations with weights that are included in the set of model
parameters $\theta$. In Section VI.1, we will describe in more detail how the
forward modeling is implemented in our present analysis framework.

Because we are analyzing event counts, it is trivial to write down the
likelihood as a product of Poisson distributions:

$\mathcal{L}_{D}(\theta)=\prod_{\text{bins}}\frac{\nu^{n}e^{-\nu}}{n!}$ (12)

where $n$ is the number of events observed in a bin and $\nu$ is that
predicted by the forward model, using a particular set of model parameter
values $\theta$.

Analyzing binned event counts has the advantage that the data points are
uncorrelated and have known probability distributions and that there are no
problems with including sparsely populated or empty bins. To our knowledge,
these features are currently unique to our approach.

##  IV Technical Developments

We now describe the main numerical and programming steps carried out while
developing our analysis code.

###  IV.1 The SPLINT package

Integrating the cross sections with standard numerical methods like two-
dimensional Gauss quadrature turns out to be quite CPU-time-consuming because
a sizable sample of NNLO structure functions has to be computed for each
integration. To solve this problem, we added to the qcdnum distribution the
splint package to construct cubic interpolation splines of structure functions
and cross sections. Sampling from splines is much faster than ab-initio
computation, while splint provides fast routines for cubic spline integration.

A cubic interpolation spline $S(u)$ in the variable $u$ is a piecewise cubic
polynomial defined on a strictly ascending set of node-points $\\{u_{i}\\}$.
The four polynomial coefficients in each node-bin are adjusted such that the
spline coincides at each node $u_{i}$ with an input value $f_{i}$, and is
continuous in the first and second derivative; the third derivative is allowed
to be discontinuous at the node points. These conditions lead to a set of
linear equations in the coefficients, which can be solved if boundary
conditions are given on the slopes at the two end points of the spline. We use
a spline algorithm that fixes the third derivatives to values estimated from
divided differences. A one-dimensional spline $S(u)$ becomes a two-dimensional
spline $S(u,v)$ by parameterizing the $u$-coefficients as cubic splines in
$v$.

The splint package constructs splines on a selected set of qcdnum $x$-$Q^{2}$
evolution grid-points. Alignment of the nodes to the evolution grid avoids
qcdnum interpolation of a user-given input function $f(x,Q^{2})$ when
constructing the spline. A coarser node-grid gives a faster spline
construction but also a less precise spline approximation so that some tuning
is necessary to balance speed versus accuracy, see Section IV.2.

Spline interpolation becomes interesting when many samples are needed of
functions that are expensive to compute, such as structure functions and cross
sections. For this purpose splint has a special routine for structure function
input that makes use of the very fast list-processing capabilities of qcdnum.
Apart from creating splines, the splint package also has routines to integrate
these over rectangular bins taking, if necessary, into account the kinematic
limit $Q^{2}\leq xs$, where $\sqrt{s}$ is the centre-of-mass energy of the
$\mathrm{e}^{\pm}\mathrm{p}$ collisions.

Inside splint, the splines are functions of the internal qcdnum variables
$u=-\ln x$ and $t=\ln Q^{2}$, which introduces a Jacobian $e^{-u}e^{t}$ in the
integrals. As a consequence, spline integration is expressed in terms of the
fundamental ($\gamma$-function) integrals

$E^{\pm}(z,n)=\int_{0}^{z}w^{n}e^{\pm w}\,\,\mathrm{d}w.$ (13)

Partial integration of Eq. 13 yields simple recursion relations between the
$E^{\pm}$ at successive values of $n$ so that they can be computed rapidly for
all powers 0–3 needed for the term-by-term integration of a cubic polynomial
inside a node-bin. When a node-bin crosses the kinematic limit, we have to
integrate both over rectangles and right-angled triangles in the $u$-$t$
plane. The triangular domains are handled by splint integration over $u$ and
Gauss integration over $t$.

To validate the procedure, we integrated splines over arbitrary rectangles,
with and without crossing the kinematic limit, both with splint and with a
two-dimensional Gauss integration routine. In all cases, the relative
difference was found to be $<10^{-9}$ with splint running about a factor of
300 faster than Gauss.

For a more detailed description of the integration algorithm in splint we
refer to the write-up, which is available from the qcdnum website.333
https://www.nikhef.nl/user/h24/qcdnum

###  IV.2 The tuning of QCDNUM and SPLINT

In the range $x>10^{-3}$ and $100<Q^{2}<3\times 10^{4}\ \mathrm{GeV}^{2}$ we
have tuned the qcdnum evolution grid and the splint spline nodes with the aim
to obtain a relative accuracy on the cross sections of better than $5\times
10^{-4}$, with the code still running at an acceptable speed. For this tuning,
the cross sections are computed at $\sqrt{s}=300\ \mathrm{GeV}$.

In Botje (2011) it is recommended to run the qcdnum evolution on a $100\times
50$ $x$-$Q^{2}$ grid with break points at $x=\\{0.2,0.4,0.6,0.8\\}$ where the
grid density doubles towards larger $x$ at each break point. Comparison with
evolution on a $300\times 150$ grid shows a relative PDF accuracy of better
than $10^{-4}$, except at very large $x$ where the parton distributions
vanish.

Because structure functions are slowly varying in $x$ and $Q^{2}$, a coarse
$22\times 7$ node grid is sufficient to yield interpolation splines with a
relative accuracy better than $5\times 10^{-4}$, and the same is true for the
cross section spline on a $100\times 25$ node grid, see Fig. 1. The latter
grid must be dense since the cross section is strongly varying.

Figure 1: Estimate of the relative error on the differential cross section
versus $x$ at $Q^{2}=100$ GeV2. Above $x=0.95$ the cross section vanishes, and
a relative error becomes ill-defined.

In Table 1, we show the CPU time needed when running the tuned code on a 2018
MacBook Pro with an Intel processor. Starting from a set of input PDFs,
computing 429 integrated cross sections takes less than 10 ms on such a
machine.

Table 1: CPU cost of computing integrated cross sections.  Subtask | Grid | CPU time [ms]
---|---|---
Evolution | $100\times 50$ | 3.6
Structure function splines (6$\times$) | $22\times 7$ | 2.9
Cross section spline | $100\times 25$ | 2.2
Integration over 429 bins |  | 0.8

###  IV.3 The QCDNUM.jl wrapper

We have developed the QCDNUM.jl package444
https://github.com/cescalara/QCDNUM.jl which is the julia interface to the
qcdnum fast QCD evolution and convolution routines. The qcdnum program is
written in fortran77 and the interface gives us access to this fast, versatile
and well-tested software in the bat analysis framework (see Section II.4).

With QCDNUM.jl, we provide a lightweight but easily usable wrapper, with all
the original qcdnum functionality available to a julia programmer. Example
programs are provided in the documentation, allowing those familiar with
qcdnum to quickly adapt their code. As part of the julia interface, we also
provide a well-documented, high-level interface to implementing and
visualizing the evolution of PDFs that may be appealing to both new and old
users of qcdnum.

The implementation in julia offers several further advantages. Thanks to
BinaryBuilder.jl 555 https://binarybuilder.org , qcdnum is automatically
compiled and installed in a cross-platform manner, with no actions needed from
the user. It is also possible to interface code relying on qcdnum with julia’s
rich functionality and ever-growing modern package ecosystem. The existing
tools for working with julia code, such as Jupyter notebooks, enable more
accessible and reproducible code. As such, the QCDNUM.jl package is a useful
addition to the community bringing the use of qcdnum to a broader audience and
modern coding practices.

###  IV.4 PartonDensity.jl package

We implement a full forward model and interface to bat in PartonDensity.jl 666
https://github.com/cescalara/PartonDensity.jl to enable the simulation of data
and inference of PDF parameters. Included are the interfaces to qcdnum and
splint, as well as the different PDF parameterizations described in this work.
We also implement the forward modeling transfer matrices from the ZEUS
experiment, described in Section VI.1. Within this framework, we provide prior
and likelihood definitions that can be used with bat, as well as tools for
visualization. The software is well documented and designed in a modular way
to allow extension to other PDF parameterizations and data sets of interest.

Several forms of standard output are available from our code. One useful
feature is a text-file summary of the posterior distribution, an example of
which is shown in the upper panel of Fig. 2.

Figure 2: Upper panel: Example of a statistical summary of a 26-dimensional
posterior distribution. Lower panel: Parameter correlation distributions for a
subset of the model parameters. The marginal distributions are shown on the
diagonal. Taken from the pseudo-data analysis described in Section VI.

The marginalized posteriors and global mode values are given for all
parameters. Many graphical representations of the output are available in
standard form, and others can be easily generated. As an example, we show in
the lower panel of Fig. 2 a two-by-two correlation plot that can be made for
any subset of the parameters.

##  V Parton Density parameterizations

The PDF determination procedure requires a well-defined set of
parameterizations for the different parton densities at the input scale
$Q^{2}_{0}$ of the evolution. A variety of forms are currently used by the
different fitting teams Cridge  _et al._ (2022); Bailey  _et al._ (2021); Hou
_et al._ (2021); Dulat  _et al._ (2016); Alekhin  _et al._ (2018, 2014); Ball
_et al._ (2022, 2017); Abt  _et al._ (2022); Abramowicz  _et al._ (2015) which
are all of the type:

$xf_{i}(x)=A_{i}x^{\lambda_{i}}F_{i}(x)(1-x)^{K_{i}}.$ (14)

Each parton species $i$ has its own set of parameters and function $F$. The
behavior of a PDF as $x\rightarrow 1$ is largely controlled by its parameter
$K$, while $\lambda$ controls the behavior as $x\rightarrow 0$. The function
$F(x)$ interpolates between these two extreme regions, and the parameter $A$
fixes the normalization.

In Section III.1 it is mentioned that the quark and gluon distributions
$xf_{i}(x)$ are required to satisfy the momentum sum rule and quark counting
rules as given in Eqs. 5 and 6, respectively. We therefore wish to construct
parameterizations that are flexible, positive-definite, easily evaluated, and
quickly integrated. We also wish, at this stage in the development of the
Bayesian framework, to restrict the number of PDF parameters as much as
possible. Given these requirements, we have fully implemented in
PartonDensity.jl a Beta-distribution parameterization that provides a suitable
starting point in a future exploration of a variety of alternatives.

###  V.1 Parameterization based on Beta distributions

To parameterize the quark densities, it is convenient to write them as valence
($q^{\mathrm{v}}$) and sea ($q^{\mathrm{s}}$) distributions,

$q+\bar{q}=(q-\bar{q})+2\bar{q}=q^{\mathrm{v}}+q^{\mathrm{s}}.$

We parameterize the valence momentum densities as

$xq^{\mathrm{v}}_{i}(x,Q^{2}_{0})=\left\\{\begin{array}[]{ll}A_{i}\;x^{\lambda_{i}}(1-x)^{K_{i}}&\
\text{for}\ i=u,d\\\ 0&\ \text{otherwise}.\end{array}\right.$ (15)

The integral of the number density, $q^{\mathrm{v}}_{i}(x,Q^{2}_{0})$ is
finite for $\lambda_{i}>0$. Similarly, for the antiquark distributions, we
have

$x\bar{q}_{i}(x,Q^{2}_{0})=A_{i}\;x^{\lambda_{\bar{q}}}(1-x)^{K_{\bar{q}}}\quad\text{for}\quad
i=\bar{u},\bar{d},\bar{s},\bar{c},\bar{b},$ (16)

where all antiquark flavors share the same $x$-dependence, but have different
normalisations $A_{i}$. For this parameterization, we require
$-1<\lambda_{\bar{q}}<0$ so that $x\bar{q}$ is integrable and increasing at
low $x$. Finally, we parameterize the gluon density as the sum of valence and
sea contributions

$\begin{split}xg(x,Q^{2}_{0})&=xg^{\mathrm{v}}(x)+xg^{\mathrm{s}}(x)\\\
&=A^{\mathrm{v}}_{g}\;x^{\lambda^{\mathrm{v}}_{g}}(1-x)^{K_{g}}+A^{\mathrm{s}}_{g}\;x^{\lambda^{\mathrm{s}}_{g}}(1-x)^{K_{\bar{q}}}.\quad\end{split}$
(17)

where we set $-1<\lambda^{\mathrm{s}}_{g}<\lambda^{\mathrm{v}}_{g}$ and set
further restrictions on these parameters in order to obtain an integrable
gluon density with a valence (sea) contribution that decreases (increases)
towards low $x$. Here $K_{g}$ is an independent parameter, while $K_{\bar{q}}$
is shared with the $x\bar{q}_{i}$ densities defined in Eq. 16.

For a Beta-distribution $xf(x)=A\,x^{\lambda}(1-x)^{K}$ we can replace the
normalization constant by the momentum integral through the relation

$\Delta=A\,\frac{\Gamma(\lambda+1)\Gamma(K+1)}{\Gamma(\lambda+K+2)}.$ (18)

The valence sum rule, introduced in Eq. 6, relates the normalizations to the
shape parameters by

$A_{i}=N^{\mathrm{v}}_{i}\;\frac{\Gamma(\lambda_{i}+K_{i}+1)}{\Gamma(\lambda_{i})\Gamma(K_{i}+1)}\qquad
i=u,d\;,$ (19)

with $N^{\mathrm{v}}_{u}=2$ and $N^{\mathrm{v}}_{d}=1$.

Using the property $\Gamma(z+1)=z\Gamma(z)$ we find from Eqs. 18 and 19 the
following relation between the total momentum $\Delta$ carried by the $u$ or
$d$ valence quarks and the shape parameters

$\Delta_{i}=N^{\mathrm{v}}_{i}\;\frac{\lambda_{i}}{\lambda_{i}+K_{i}+1}\qquad
i=u,d.$ (20)

###  V.2 Prior parameter constraints

Given a set of free PDF parameters, the challenge is to include as many
physically motivated constraints on the priors as possible.

For this, it is advantageous to replace the normalization constants, $A$,
which are not straightforward for interpretation, with the momentum fractions,
$\Delta$, either by numerically integrating the PDFs or by using Eq. 18 in
case the PDFs are parameterized in terms of Beta-distributions.

In that case Eq. 20 offers two alternatives for the parameters of the up and
down valence distributions: (i) leave the shape parameters $\lambda$ and $K$
free and thereby fix the momentum fraction $\Delta$ or (ii) leave $\Delta$
free and thereby fix one of the two shape parameters.

The first alternative allows to specify the shape of the valence densities,
but it is not trivial to ensure that $\Delta_{d}+\Delta_{u}<1$ and
$\Delta_{d}<\Delta_{u}$ over the full support of the prior.

In the analysis of the ZEUS high-$x$ data, we have chosen to fix the low-$x$
shape parameter $\lambda$ and leave free the high-$x$ shape parameter $K$ and
the momentum fraction $\Delta$. In this way we include in the model parameters
the complete set of PDF momenta, subject to the sum rule constraint of Eq. 5
given in Section III.1.

It is convenient to use a Dirichlet distribution as a joint prior for the
momenta Betancourt (2012). A Dirichlet distribution ${\rm
Dir}(\alpha_{1},\ldots,\alpha_{k})$ of $k$ independent variables
$u_{i}\in[0,1]$ lives on a $(k-1)$-dimensional manifold defined by $\sum
u_{i}=1$ so that the momentum sum rule is automatically satisfied. It is a
multivariate generalization of the Beta distribution; for instance ${\rm
Beta}(\alpha_{1},\alpha_{2})$ of one variable $u$ is the same as ${\rm
Dir}(\alpha_{1},\alpha_{2})$ of two variables $(u_{1},u_{2})$ with
$u_{1}+u_{2}=1$.

With a suitable choice of the shape parameters $\alpha$, it is possible to
satisfy expectations such as that, asymptotically, gluons and quarks carry
approximately the same momentum, that valence up-quarks carry about twice the
momentum of valence down-quarks, and that the heavier quarks carry little
momentum.

The spectator counting rules of Brodsky and Farrar Brodsky and Farrar (1973)
give a first expectation for the ranges of the shape parameters
$K_{u},K_{d},K_{\bar{q}}$ and $K_{g}$. Furthermore, a body of PDF results
available from the literature indicate the preferred range for these
parameters.

##  VI Validation

To validate our Bayesian tools, we analyzed sets of simulated data and
compared the posterior distributions to the parameter input values. Because
the tools were initially developed for the analysis Aggarwal  _et al._ (2023)
of the ZEUS high-$x$ and large-$Q^{2}$ data Abramowicz  _et al._ (2014), our
simulations were restricted to the kinematic range covered by these data. An
adequate representation is obtained by using the simple Beta parameterization
described in the previous sections, at an input scale of $Q^{2}_{0}=100$ GeV2.
Some details of the ZEUS data are given below.

###  VI.1 ZEUS data and systematics

In Abramowicz  _et al._ (2014), the ZEUS collaboration has published
$\mathrm{e}^{\pm}\mathrm{p}$ deep-inelastic scattering data covering the range
$0.03\leq x\leq 1$ and $650\leq Q^{2}\leq 20000$ GeV2. These data are unique
in providing measurements up to $x=1$ in the high $Q^{2}$ regime where higher-
twist effects are absent so that a clean analysis can be performed based on
the perturbative QCD evolution equations. In Aggarwal  _et al._ (2023), it is
shown that the data mainly constrain the parameters of the valence up-quark,
as will also become apparent in the figures and results given below.

The data are presented as counts in 153 bins, separately for the e-p and e+p
data sets.

To enable forward modeling, ZEUS has provided in Abt  _et al._ (2020) the
transfer matrix $A$ as defined by Eq. 11 in Section III.2. This matrix is
written as the product of two matrices, $R$ and $T$, which account for
radiative and reconstruction effects, respectively. Also provided are
variations on $T$ to enable the evaluation of systematic uncertainties. The
QED-Born level binning of 429 bins is finer than the experimental binning.

The $R$ matrices are $429\times 429$ diagonal and correct the Born-level cross
sections for $O(\alpha)$ QED effects.

A reconstruction-matrix element $T_{ij}$ gives the probability that an event
generated in $x$-$Q^{2}$ bin $j$ is reconstructed in experimental bin $i$ so
that $T$ has a dimension of $153\times 429$.

The uncertainties that need to be accounted for in the ZEUS data are the
following.

  1. 1.

The uncertainty in the luminosity of each of the $\mathrm{e}^{+}\mathrm{p}$
and $\mathrm{e}^{-}\mathrm{p}$ data sets;

  2. 2.

The uncertainties related to the imperfect understanding of detector effects
such as acceptance, energy resolution, etc. These are completely correlated
between the two data sets. Variational matrices $T^{\prime}$ are provided for
8 sources of uncertainty in each of the $\mathrm{e}^{-}\mathrm{p}$ and
$\mathrm{e}^{+}\mathrm{p}$ data sets.

It turns out that the ZEUS luminosity uncertainty of 1.8% is the most
important source of uncertainty in this analysis.

From the matrices provided, we can construct for each source of systematic
error the deviation matrix

$A^{\prime}=R\,(T^{\prime}-T)$

which corresponds to a one standard deviation uncertainty, assuming that the
systematic errors are Gaussian-distributed. Denoting the one standard
deviation uncertainty on the luminosity by $\beta_{0}$, and including weight
factors $\beta_{0},\ldots,\beta_{8}$ in the set of model parameters $\theta$,
Eq. 11 can be written as, in matrix notation,

$\nu=\mathcal{L}\,(1+0.018\cdot\beta_{0})\left(A+\sum_{k}\beta_{k}A^{\prime}_{k}\right)\lambda,$
(21)

where the sum runs over the 8 sources of systematic error and the $\beta$
parameters are taken to be normally distributed with unit width.

###  VI.2 Model data

To validate the analysis framework, we generated pseudo-data by computing
event predictions $\nu$ and event counts $n$ that are Poisson distributed with
mean $\nu$. With the input parameters given in Table 2

Table 2: Parameter values used in the data simulation.

$\bm{\Delta}\times 10^{3}$
---
$u^{V}$ | $d^{V}$ | $g^{V}$ | $g^{S}$ | $2\bar{u}$ | $2\bar{d}$ | $2\bar{s}$ | $2\bar{c}$ | $2\bar{b}$
228 | 104 | 249 | 249 | 104 | 52 | 10 | 5 | 0.5

$K_{u}$ | $K_{d}$ | $\lambda^{\mathrm{v}}_{g}$ | $\lambda^{\mathrm{s}}_{g}$ | $K_{g}$ | $\lambda_{\bar{q}}$ | $K_{\bar{q}}$ |  $\bm{\beta}$
---|---|---|---|---|---|---|---
3.70 | 3.70 | 0.50 |  $-$0.50 | 5.0 |  $-$0.50 | 6.0 | 0

the simulation fairly well describes the actual ZEUS data. All $\beta$
parameters are set to zero so that there are no systematic biases in the model
data. To study convergence, we have also generated sets with the ZEUS
luminosity increased by a factor of 50 so that the likelihood becomes almost a
delta function.

###  VI.3 Fits to the model data

Unless otherwise stated, we use the priors listed in Table 3 in the fits777 In
a Bayesian analysis, parameters are not fitted in the usual sense. In the
context of this paper, a fitted (or free) parameter is one that has a prior
_distribution_ assigned to it. Fixed parameters, on the other hand, have
single-valued priors, thereby reducing the dimension of the parameter space to
be explored.  to the pseudo-data.

Table 3: Priors used in the analysis of the pseudo-data sets. There are 9 parameters in the vector $\bm{\Delta}$ and 10 in $\bm{\beta}$. The normal distributions are truncated to the range indicated, and their mean and standard deviation are given in brackets. | Prior | Range
---|---|---
$\bm{\Delta}$ | Dir(20, 10, 20, 20, 5, 2.5, 1.5, 1.5, 0.5) | $[0,1]$
$K_{u}$ | Normal(3.5, 0.5) | $[1,6.5]$
$K_{d}$ | Normal(3.5, 0.5) | $[1,6.5]$
$\lambda^{\mathrm{v}}_{g}$ | Uniform | $[0,1]$
$\lambda^{\mathrm{s}}_{g}$ | Uniform | $[-1,-0.1]$
$K_{g}$ | Normal(4, 1.5) | $[1,8.5]$
$\lambda_{\bar{q}}$ | Uniform | $[-1,-0.1]$
$K_{\bar{q}}$ | Normal(4, 1.5) | $[1,9.5]$
$\bm{\beta}$ | Normal(0, 1) | $[-5,5]$

The matrices for describing the transformation from QED Born-level cross
sections to the observed level are taken from the ZEUS Collaboration Abt  _et
al._ (2020) for all tests described in this paper. Tests were performed both
with the systematic $\beta$ parameters fixed to zero, or with them left free.

Fig. 3

Figure 3: The prior and posterior probability distributions of
$(\Delta_{u},K_{u})$ (upper) and $(\Delta_{d},K_{d})$ (lower) from the nominal
(left) and high-luminosity (right) pseudo-data sets. The red crosses indicate
the true parameter values.

shows two-dimensional prior and posterior distributions for the parameters
$(\Delta_{u},K_{u})$ and $(\Delta_{d},K_{d})$ obtained from the nominal and
high-luminosity pseudo-data sets. As we have already observed in Aggarwal  _et
al._ (2023), a comparison of posterior to prior shows a very significant
knowledge update for the valence up-quark parameters, and less so for those of
the valence down-quark. The knowledge update on the parameters of the gluon
and sea distributions (not shown) is also rather limited, as is expected from
an analysis of high-$x$ data.

Comparison of the left- and right-hand plots in Fig. 3 shows that the
posterior converges to the true values, indicating that there is no
appreciable bias introduced by our analysis procedure. We have verified that
this is the case for all fitted parameters, although convergence with
increasing statistics is slow or absent for those that are ill-constrained by
the data.

Parameters that are well-constrained by the data should be relatively
insensitive to the choice of prior. Many different data sets were analyzed
using shifted prior choices to determine if these significantly affected the
results. Two such tests for the up-valence parameters are shown in Fig. 4.

Figure 4: The posterior probability distributions of $(\Delta_{u},K_{u})$ from
two analyses with strongly biased priors (see text). True values are indicated
by the red crosses. The legend is the same as that of Fig. 4.

In the left-hand plot, the momentum fraction $\Delta_{u}$ was strongly biased
upward borrowing momentum from the gluon while the prior for $K_{u}$ also had
an upward bias. In the right-hand plot, the gluon momentum was strongly biased
upward, together with a downward bias of $K_{u}$. In both cases, the posterior
well reproduced the true values of $\Delta_{u}$ and $K_{u}$.

As a further example of the information made available by the full posterior
probability density we plot in Fig. 5

Figure 5: Correlations between parton momenta from the analysis of the high-
luminosity pseudo-data. True values are indicated by the red crosses. The
legend is the same as that of Fig. 4.

various correlations among the momentum components of the proton obtained from
the analysis of the high-luminosity data set. Here $\Delta_{\rm sea}$ is the
fractional momentum of the sea (anti-)quarks, summed over all quark flavors.
Again, it is seen that the true values are nicely reproduced and also that the
momentum carried by the valence up quark is very well constrained and weakly
sensitive to the other components. The momenta carried by the sea quarks and
the gluon density are anti-correlated, as is to be expected since the sea is
generated from gluon splitting. For another correlation plot, we refer to Fig.
2 in Section IV.4.

The results shown up to now were obtained from fits where the systematic
parameters were kept fixed to zero. Leaving these parameters free showed that
the data hardly constrain them and that they remain un-correlated, except for
the $\mathrm{e}^{+}\mathrm{p}$ and $\mathrm{e}^{-}\mathrm{p}$ luminosity
parameters that show a strong correlation, as expected. However, a fit with
free systematic parameters is always preferable since that allows for error
propagation by marginalizing the posterior over them. However, with the data
sets studied here, this had a minor effect.

We mentioned in Section V.2 an alternative scheme where, instead of fitting
all parton momenta, those of the up and down valence quarks are fixed by
fitting, instead, their low-$x$ $\lambda$ shape parameters. In such a fit, we
have removed the priors on $\Delta_{u,d}$ and introduced those on
$\lambda_{u,d}$ as normal distributions of unit width centered at $0.5$ and
truncated to a range $[0.2,0.8]$. It turned out that the results were very
similar to those from the standard parameterization.

Figure 6: True parton momentum densities (red curves) compared to those
computed from 100 samples of the posterior distribution (gray curves) from
data at nominal luminosity.

In Fig. 6 (note the vertical logarithmic scale) we show input parton
distributions used in the generation of the pseudo-data compared to those
computed from 100 parameter sets randomly sampled from the posterior
distribution of the data set at nominal luminosity. Again, it is seen that the
valence up-quark distribution is faithfully reproduced. The valence down-quark
distribution is poorly constrained and has a wide spread. The general features
of the gluon and up-antiquark distributions are correctly reproduced.

##  VII Goodness-of-Fit Tests

A standard Bayesian analysis does not provide goodness-of-fit tests that
compare a single model to data. In fact, a single-model Bayesian analysis
yields, as a result, the posterior probability distribution of the model
parameters, and no statement is made concerning the validity of the model
itself. To investigate model dependence, a choice is to be based on posterior
probabilities calculated for several alternative models F. Beaujean (2011).
Such model selection is beyond the scope of this paper and will be used in a
future extension of the framework to investigate the sensitivity to different
parameterizations of the input PDFs.

Although not strictly Bayesian, it is possible, and maybe welcome, to provide
a goodness-of-fit test based on a single posterior probability distribution.
We consider two possibilities below.

###  VII.1 Posterior predictive check

In a posterior predictive check, the actual data are compared to pseudo-data
distributions generated from parameter sets that are sampled from the
posterior distribution. The quality of the model can then be judged by
observing good or bad overlap of the pseudo- and observed data.

A graphical example of such a test is shown in Fig. 7

Figure 7: Event number distributions computed with parameters sampled from the
posterior distribution (bands) compared to the simulated
$\mathrm{e}^{-}\mathrm{p}$ scattering data (dots). For clarity, the box at
large $x$ is shown enlarged in the inset.

which shows good agreement. This is not surprising since the data themselves
were generated from a given parameter set, as is described in Section VI. Here
the figure confirms that the input parameters are well-reproduced and that the
analysis framework does not introduce significant bias in the result.

###  VII.2 Posterior mode $\chi^{2}$ test

A simple goodness-of-fit test is provided by computing the Pearson $\chi^{2}$,
defined by

$\chi^{2}_{P}=\sum_{i}\frac{(n_{i}-\nu_{i})^{2}}{\nu_{i}},$

where the sum runs over all bins $i$ with observed event numbers $n_{i}$, and
expected number of events $\nu_{i}$ as calculated from the posterior global
mode parameter values. In the presence of a sizable amount of sparsely
populated or empty bins, $\chi^{2}_{P}$ does not follow a standard $\chi^{2}$
distribution. In this case, it is better to calculate $\chi^{2}_{P}$ from a
large set of simulated pseudo-data. In Fig. 8

Figure 8: The $\chi^{2}$ distribution computed from the combined
$\mathrm{e}^{+}\mathrm{p}$ and $\mathrm{e}^{-}\mathrm{p}$ analysis results.
The curve shows the prediction for 306 degrees of freedom.

we show a histogram of the $\chi^{2}_{P}$ distribution obtained from many toy-
model $\mathrm{e}^{\pm}\mathrm{p}$ data sets that have $\nu_{i}$ fixed and
$n_{i}$ Poisson distributed around $\nu$. From such histograms, it is
straightforward to compute $p$-values by normalizing the histogram and summing
the bin contents above the observed value of $\chi^{2}_{P}$. In our tests, we
have found a fairly flat distribution of $p$-values, as expected.

We have seen in our studies that the maximum posterior probability results do
not coincide with a minimum value of $\chi^{2}_{P}$. This is not surprising
since a maximum posterior is not necessarily a maximum likelihood.
Nevertheless, the distribution in Fig. 8 closely resembles that of a standard
$\chi^{2}$ and, as explained in F. Beaujean (2011); Caldwell (2019), small
$p$-values can be understood from a Bayesian perspective as implying that the
model is open to improvement.

##  VIII Summary and Outlook

We have developed a novel parton density analysis code that allows for a full
Bayesian posterior probability determination and also supports a forward
modeling approach. The open-source code has been thoroughly tested and is now
available for distribution. In this paper, we report on its structure, the
technical developments that have been made in realizing the code, and a series
of validation tests that have been performed. We believe that the code is
reliable, well-documented, and easy to use.

To date, the code has been used exclusively for the analysis of high-$x$ and
high-$Q^{2}$ $e^{\pm}p$ deep inelastic scattering data. We look forward to
also extending the analysis to other data sets, including those reported as
differential cross sections at the QED Born level, although such an analysis
cannot benefit from the statistical rigor offered by the forward modeling
approach. In fact, the information needed to enable such an approach is often
not made available in proper form by the experiments.

An important step in achieving this is to make the analysis code run much
faster. Here there is ample room for improvement by parallelizing computations
through threading or forking, by improving the MCMC sampling efficiency and by
speeding up the QCD evolution of the PDFs.

We also intend to extend the framework to investigate more flexible PDF
parameterizations using Bayesian model selection techniques.

## Acknowledgements

The authors would like to thank Richard Hildebrandt for extensive
investigations of alternative PDF parameterizations, which we intend to report
in a future publication.

R. Aggarwal acknowledges the support of Savitribai Phule Pune University. F.
Capel was employed by the Excellence Cluster ORIGINS, which is funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094-390783311, during much of the project
duration.

## References

  * Cridge  _et al._ (2022) T. Cridge _et al._ , Eur. Phys. J. C 82, 90 (2022), arXiv:2111.05357 [hep-ph] .
  * Bailey  _et al._ (2021) S. Bailey _et al._ , Eur. Phys. J. C 81, 341 (2021), arXiv:2012.04684 [hep-ph] .
  * Hou  _et al._ (2021) T.-J. Hou _et al._ , Phys. Rev. D 103, 014013 (2021), arXiv:1912.10053 [hep-ph] .
  * Dulat  _et al._ (2016) S. Dulat _et al._ , Phys. Rev. D 93, 033006 (2016), arXiv:1506.07443 [hep-ph] .
  * Alekhin  _et al._ (2018) S. Alekhin, J. Blümlein, and S. Moch, Eur. Phys. J. C 78, 477 (2018), arXiv:1803.07537 [hep-ph] .
  * Alekhin  _et al._ (2014) S. Alekhin, J. Blumlein, and S. Moch, Phys. Rev. D 89, 054028 (2014), arXiv:1310.3059 [hep-ph] .
  * Ball  _et al._ (2022) R. D. Ball _et al._ (NNPDF), Eur. Phys. J. C 82, 428 (2022), arXiv:2109.02653 [hep-ph] .
  * Ball  _et al._ (2017) R. D. Ball _et al._ (NNPDF), Eur. Phys. J. C 77, 663 (2017), arXiv:1706.00428 [hep-ph] .
  * Abt  _et al._ (2022) I. Abt _et al._ (H1, ZEUS), Eur. Phys. J. C 82, 243 (2022), arXiv:2112.01120 [hep-ex] .
  * Abramowicz  _et al._ (2015) H. Abramowicz _et al._ (H1, ZEUS), Eur. Phys. J. C 75, 580 (2015), arXiv:1506.06042 [hep-ex] .
  * Cocuzza  _et al._ (2021) C. Cocuzza _et al._ (Jefferson Lab Angular Momentum (JAM)), Phys. Rev. D 104, 074031 (2021), arXiv:2109.00677 [hep-ph] .
  * Hunt-Smith  _et al._ (2022) N. T. Hunt-Smith _et al._ , Phys. Rev. D 106, 036003 (2022), arXiv:2206.10782 [hep-ph] .
  * Bezanson  _et al._ (2017) J. Bezanson _et al._ , SIAM review 59, 65 (2017).
  * Schulz  _et al._ (2021) O. Schulz _et al._ , SN Computer Science 2, 210 (2021).
  * Aggarwal  _et al._ (2023) R. Aggarwal _et al._ , Phys. Rev. Lett. 130, 141901 (2023), arXiv:2209.06571 [hep-ph] .
  * Abramowicz  _et al._ (2014) H. Abramowicz _et al._ (ZEUS), Phys. Rev. D 89, 072007 (2014), [Erratum: Phys.Rev.D 106, 079902 (2022)], arXiv:1312.4438 [hep-ex] .
  * Abt  _et al._ (2020) I. Abt _et al._ (ZEUS), Phys. Rev. D 101, 112009 (2020), [Erratum: Phys.Rev.D 106, 079901 (2022)], arXiv:2003.08742 [hep-ex] .
  * Betancourt (2012) M. Betancourt, AIP Conference Proceedings 1443, 157 (2012).
  * Gribov and Lipatov (1972a) V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 438 (1972a).
  * Gribov and Lipatov (1972b) V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 675 (1972b).
  * Lipatov (1975) L. N. Lipatov, Sov. J. Nucl. Phys. 20, 94 (1975).
  * Dokshitzer (1977) Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977).
  * Altarelli and Parisi (1977) G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977).
  * Botje (2011) M. Botje, Comput. Phys. Commun. 182, 490 (2011), arXiv:1005.1481 [hep-ph] .
  * Tarasov  _et al._ (1980) O. V. Tarasov, A. A. Vladimirov, and A. Y. Zharkov, Phys. Lett. B 93, 429 (1980).
  * Curci  _et al._ (1980) G. Curci, W. Furmanski, and R. Petronzio, Nucl. Phys. B 175, 27 (1980).
  * Furmanski and Petronzio (1980) W. Furmanski and R. Petronzio, Phys. Lett. B 97, 437 (1980).
  * Furmanski and Petronzio (1982) W. Furmanski and R. Petronzio, Z. Phys. C 11, 293 (1982).
  * Larin and Vermaseren (1993) S. A. Larin and J. A. M. Vermaseren, Phys. Lett. B 303, 334 (1993), arXiv:hep-ph/9302208 .
  * Chetyrkin  _et al._ (1997) K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, Phys. Rev. Lett. 79, 2184 (1997), arXiv:hep-ph/9706430 .
  * Buza  _et al._ (1998) M. Buza _et al._ , Eur. Phys. J. C 1, 301 (1998), arXiv:hep-ph/9612398 .
  * Moch  _et al._ (2004) S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B 688, 101 (2004), arXiv:hep-ph/0403192 .
  * Vogt  _et al._ (2004) A. Vogt, S. Moch, and J. A. M. Vermaseren, Nucl. Phys. B 691, 129 (2004), arXiv:hep-ph/0404111 .
  * Ball  _et al._ (2016) R. Ball _et al._ , Phys. Lett. B 754, 49 (2016), arXiv:1510.00009 [hep-ph] .
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Sanchez Guillen  _et al._ (1991) J. Sanchez Guillen _et al._ , Nucl. Phys. B 353, 337 (1991).
  * van Neerven and Zijlstra (1991) W. L. van Neerven and E. B. Zijlstra, Phys. Lett. B 272, 127 (1991).
  * Zijlstra and van Neerven (1991) E. B. Zijlstra and W. L. van Neerven, Phys. Lett. B 273, 476 (1991).
  * Zijlstra and van Neerven (1992) E. B. Zijlstra and W. L. van Neerven, Phys. Lett. B 297, 377 (1992).
  * van Neerven and Vogt (2000a) W. L. van Neerven and A. Vogt, Nucl. Phys. B 568, 263 (2000a), arXiv:hep-ph/9907472 .
  * van Neerven and Vogt (2000b) W. L. van Neerven and A. Vogt, Nucl. Phys. B 588, 345 (2000b), arXiv:hep-ph/0006154 .
  * Brodsky and Farrar (1973) S. J. Brodsky and G. R. Farrar, Phys. Rev. Lett. 31, 1153 (1973).
  * F. Beaujean (2011) F. Beaujean, Phys. Rev. D 83, 012004 (2011), arXiv:1011.1674 [physics.data-an] .
  * Caldwell (2019) A. Caldwell, Ann. d. Physik 531, 1700457 (2019).