* [39] Eugene Gorsky, Graham Hawkes, Anne Schilling, and Julianne Rainbolt.  Generalized $q,t$-Catalan numbers.  Algebraic Combinatorics, 3(4):855–886, 2020.
  * [40] Eugene Gorsky, Matthew Hogancamp, Anton Mellit, and Keita Nakagane.  Serre duality for Khovanov-Rozansky homology.  Selecta Math. (N.S.), 25(5):Paper No. 79, 33, 2019.
  * [41] Eugene Gorsky, Andrei Neguţ, and Jacob Rasmussen.  Flag Hilbert schemes, colored projectors and Khovanov-Rozansky homology.  Advances in Mathematics, 378:107542, 2021.
  * [42] Mikhail Gorsky.  Geometry and combinatorics of subword complexes and their dual polytopes.  Ph.D. thesis (in Russian), Steklov Mathematical Institute, 2014.
  * [43] Mikhail A Gorsky.  Subword complexes and 2-truncated cubes.  Russian Mathematical Surveys, 69(3):572–574, 2014.
  * [44] Mikhail A Gorsky.  Subword complexes and edge subdivisions.  Proceedings of the Steklov Institute of Mathematics, 286(1):114–127, 2014.
  * [45] Mark Gross, Paul Hacking, Sean Keel, and Maxim Kontsevich.  Canonical bases for cluster algebras.  Journal of the American Mathematical Society, 31(2):497–608, 2018\.
  * [46] Tamás Hausel and Fernando Rodriguez-Villegas.  Mixed Hodge polynomials of character varieties.  Inventiones mathematicae, 174(3):555–624, 2008.
  * [47] James Hughes.  Weave Realizability for D–type.  arXiv preprint arXiv:2101.10306, 2021.
  * [48] Grace Ingermanson.  Cluster algerbras of open Richardson varieties.  Ph.D. Thesis, University of Michigan, 2019.
  * [49] Tamás Kálmán.  Braid-positive Legendrian links.  Int. Math. Res. Not., pages Art ID 14874, 29, 2006.
  * [50] Rachel Karpman.  Bridge graphs and Deodhar parametrizations for positroid varieties.  Journal of Combinatorial Theory, Series A, 142:113–146, 2016.
  * [51] Mikhail Khovanov.  Triply-graded link homology and Hochschild homology of Soergel bimodules.  Internat. J. Math., 18(8):869–885, 2007.
  * [52] Mikhail Khovanov and Lev Rozansky.  Matrix factorizations and link homology. II.  Geom. Topol., 12(3):1387–1425, 2008.
  * [53] Allen Knutson, Thomas Lam, and David E Speyer.  Positroid varieties: juggling and geometry.  Compositio Mathematica, 149(10):1710–1752, 2013.
  * [54] Allen Knutson, Thomas Lam, and David E Speyer.  Projections of Richardson varieties.  Journal für die reine und angewandte Mathematik (Crelles Journal), 2014(687):133–157, 2014.
  * [55] Allen Knutson and Ezra Miller.  Subword complexes in Coxeter groups.  Advances in Mathematics, 184(1):161–176, 2004.
  * [56] Thomas Lam and David E Speyer.  Cohomology of cluster varieties. I. Locally acyclic case.  arXiv preprint arXiv:1604.06843, 2016.
  * [57] Bernard Leclerc.  Cluster structures on strata of flag varieties.  Advances in Mathematics, 300:190–228, 2016.
  * [58] Thomas H Lenagan and Milen T Yakimov.  Prime factors of quantum Schubert cell algebras and clusters for quantum Richardson varieties.  Journal für die reine und angewandte Mathematik, 2019(750):123–156, 2019.
  * [59] George Lusztig.  Total positivity in reductive groups.  In Lie theory and geometry, pages 531–568. Springer, 1994.
  * [60] George Lusztig.  Total positivity in partial flag manifolds.  Representation Theory of the American Mathematical Society, 2(3):70–78, 1998.
  * [61] Bethany Marsh and Konstanze Rietsch.  Parametrizations of flag varieties.  Representation Theory of the American Mathematical Society, 8(9):212–242, 2004.
  * [62] Anton Mellit.  Cell decompositions of character varieties.  arXiv preprint arXiv:1905.10685, 2019.
  * [63] Etienne Ménard.  Algèbres amassées associées aux variétés de Richardson ouvertes : un algorithme de calcul de graines initiales.  Theses, ComUE Normandie Université, March 2021.
  * [64] Matthew R Mills.  Maximal green sequences for quivers of finite mutation type.  Advances in Mathematics, 319:182–210, 2017.
  * [65] Greg Muller.  Locally acyclic cluster algebras.  Advances in Mathematics, 233(1):207–247, 2013.
  * [66] Greg Muller and David Speyer.  Cluster algebras of Grassmannians are locally acyclic.  Proceedings of the American Mathematical Society, 144(8):3267–3281, 2016.
  * [67] Greg Muller and David E Speyer.  The twist for positroid varieties.  Proceedings of the London Mathematical Society, 115(5):1014–1071, 2017.
  * [68] Yu Pan.  Exact Lagrangian fillings of Legendrian $(2,n)$ torus links.  Pacific J. Math., 289(2):417–441, 2017.
  * [69] Vincent Pilaud and Christian Stump.  EL-labelings and canonical spanning trees for subword complexes.  In Discrete geometry and optimization, pages 213–248. Springer, 2013.
  * [70] Vincent Pilaud and Christian Stump.  Brick polytopes of spherical subword complexes and generalized associahedra.  Advances in Mathematics, 276:1–61, 2015.
  * [71] Alexander Postnikov.  Total positivity, Grassmannians, and networks.  http://math.mit.edu/ apost/papers/tpgrass.pdf, 2006.
  * [72] Fan Qin.  Bases for upper cluster algebras and tropical points.  arXiv preprint arXiv:1902.09507, 2019.
  * [73] Konstanze Christina Rietsch.  Total positivity and real flag varieties.  PhD thesis, Massachusetts Institute of Technology, 1998.
  * [74] Joshua S Scott.  Grassmannians and cluster algebras.  Proceedings of the London Mathematical Society, 92(2):345–380, 2006\.
  * [75] Khrystyna Serhiyenko, Melissa Sherman-Bennett, and Lauren Williams.  Cluster structures in Schubert varieties in the Grassmannian.  Proceedings of the London Mathematical Society, 119(6):1694–1744, 2019.
  * [76] Linhui Shen and Daping Weng.  Cluster structures on double Bott-Samelson cells.  arXiv preprint arXiv:1904.07992, 2019.
  * [77] Vivek Shende, David Treumann, Harold Williams, and Eric Zaslow.  Cluster varieties from Legendrian knots.  Duke Math. J., 168(15):2801–2871, 2019.
  * [78] Minh-Tam Trinh.  Algebraic Braids and Geometric Representation Theory.  ProQuest LLC, Ann Arbor, MI, 2020.  Thesis (Ph.D.)–The University of Chicago.
  * [79] Ben Webster and Geordie Williamson.  A geometric model for Hochschild homology of Soergel bimodules.  Geometry & Topology, 12(2):1243–1263, 2008.


[1]Octavious Smiley

[1]Biostatistics, Harvard University, 677 Huntington Ave, Boston, 02115, MA,
USA

# Approximate Inference for Longitudinal Mechanistic HIV Contact Networks

<EMAIL_ADDRESS>Till Hoffmann<EMAIL_ADDRESS>Jukka-
Pekka Onnela<EMAIL_ADDRESS>*

###### Abstract

Network models are increasingly used to study infectious disease spread.
Exponential Random Graph models have a history in this area, with scalable
inference methods now available. An alternative approach uses mechanistic
network models. Mechanistic network models directly capture individual
behaviors, making them suitable for studying sexually transmitted diseases.
Combining mechanistic models with Approximate Bayesian Computation allows
flexible modeling using domain-specific interaction rules among agents,
avoiding network model oversimplifications. These models are ideal for
longitudinal settings as they explicitly incorporate network evolution over
time. We implemented a discrete-time version of a previously published
continuous-time model of evolving contact networks for men who have sex with
men (MSM) and proposed an ABC-based approximate inference scheme for it. As
expected, we found that a two-wave longitudinal study design improves the
accuracy of inference compared to a cross-sectional design. However, the gains
in precision in collecting data twice, up to 18%, depend on the spacing of the
two waves and are sensitive to the choice of summary statistics. In addition
to methodological developments, our results inform the design of future
longitudinal network studies in sexually transmitted diseases, specifically in
terms of what data to collect from participants and when to do so.

###### keywords:

mechanistic model, networks, HIV, ABC, inference, MSM, agent based modeling

##  1 Introduction

Networks are used to study a range of systems with interactions or
dependencies among their agents, such as the behavior in supply chains and the
stock market [1], protein-protein interactions in biological systems [2], and
disease transmission on local and global scales [3]. In the study of disease
transmission dynamics, the contact structure of a population can be naturally
represented as a network, and this representation is especially useful if the
contacts persist over time, as is often the case for sexual interactions.
Disease dynamics are then driven by interactions (represented by edges) among
susceptible and infectious individuals (represented by nodes). More generally,
many of these systems arise from stochastic processes forming or dissolving
interactions over time that must be accounted for when doing inference.

There are (at least) two main paradigms of networks models: statistical and
mechanistic. Statistical network models prioritize tractable likelihoods to
facilitate inference at the expense of model flexibility. For example, the
Erdős–Rényi graph, also known as the Bernoulli random graph, assumes that each
node pair is connected independently and with identical probability. Hence,
the likelihood of the number of edges is the standard binomial likelihood with
a fixed number of nodes and inference readily follows because a graph is
completely identified by its node and edge sets. It also follows that
Erdős–Rényi has a binomial (approximately Poisson) degree distribution.

This generative mechanism however clearly does not map well to most real-world
networks. This is easily seen in the World Wide Web (WWW). In this scenario,
each website is represented by a node and a directed connection (hyperlink)
between websites occurs when one links to the other. Unlike Erdős–Rényi
networks, where the degree distribution follows a binomial distribution, the
degree distribution here follows a power-law where more successful websites
tend to grow their connections faster than others [4]. Exponential Random
Graph Models (ERGMs) are generalizations of the Erdős–Rényi model. They
represent a probability distribution of graphs on a fixed node set, where the
probability of observing a graph is dependent on the presence of the various
configurations specified by the model [5]. A typical graph in this
distribution can be interpreted as the aggregate of the local configurations,
and slight errors in estimating the local configuration counts can alter
beliefs about the distribution [6].

Mechanistic models assume that the observed network is generated by a small
set of mechanistic rules. The canonical example is the Barabási-Albert (BA)
model. Nodes are added one by one to a growing network and each node connects
to $m$ previously existing nodes with probability proportional to the nodes’
current degree [7]. This so-called preferential attachment mechanism readily
generates power-law degree distributions, which are a type of broad-tailed
degree distribution that are characteristic of many empirical networks,
including that of the WWW. Apart from the target number of nodes, $n$, the
classic BA model only has one free parameter, $m$. In this case, the fully
grown graph has approximately $nm$ edges, and $m$ can be inferred by dividing
the number of edges by the number of nodes $n$, i.e., $m$ is approximately
equal to the average degree of the network. However, even for moderately
complex models, the likelihood of the full network becomes intractable due to
the fact that the insertion order of the nodes is (usually) not known. Because
the graph is sequentially dependent on the previous iteration as it grows, the
number of possible graph realizations grows exponentially with the number of
added nodes.

Networks have provided insights to major public health problems such as the
spread of HIV [8], the opioid crisis [9], and interventions with people who
inject drugs [10]. Wertheim et al. noted that HIV is an evolving disease and
constructed a disease transmission network using gene sequencing by tracking
the evolutionary path of the virus and inferring edges by measuring the
similarity of the virus within different individuals. Using an inferred
transmission network, they developed a network statistic that was able to
detect community level effects of HIV in a clinical trial setting that could
help thwart future infections [8]. Aroke et al. showed the benefit of peer
influence and concluded that individuals who have a diagnosis of opioid use
disorder or use many prescribers may help promote positive health behaviors in
an opioid prescription network due to the influence of their direct peers on
the network structure. They came to this conclusion by showing the type of
opioid that an individual uses and their number of prescribers were identified
as significant predictors of high betweenness centrality giving them influence
over the network at large [9]. Rolls et al. model network data involving
people who inject drugs, using validation techniques, so that these networks
can be simulated and intervention strategies could be explored [10].

There are several mechanistic models for studying the impact of men who have
sex with men (MSM) contact networks and their impact on HIV transmission [11,
12, 13]. Birkett et al. used a data-driven simulation model to understand the
impact of network-level mechanisms and STI infections on the spread of HIV
among Young Men who have Sex with Men (YMSM) [11]. Mei et al. introduced the
concept of a Complex Agent Network (CAN) to model the HIV epidemics by
combining agent-based modelling and complex networks [12]. An especially
interesting model was introduced by Hansson et al. to study the role of casual
contacts on the HIV epidemic in Stockholm, Sweden. Their research was used to
recommend interventions to reduce transmission rates [13]. Padeniya et al.
notes Hansson and others’ contribution to intervention strategies as they
sought to mathematically model the role of female-sex-worker-client
interactions for gonorrhoea transmission [14]. Vajdi et al. noted Hansson’s
choice to model instantaneous casual relationships, and investigated a dynamic
model for casual relationships, a two-layer temporal network model, and SIS
mean-field equations [15]. A common approach for inference in these papers is
to propose mechanisms for contact formation, simulate the spread of disease on
the network, and modify parameter values to match disease prevalence to that
observed in their respective populations without directly validating their
mechanism.

Most scientific studies involving human subjects can be divided into cross-
sectional and longitudinal. In cross-sectional studies, measurements are
obtained at only a single point in time. The distinguishing feature of
longitudinal studies is that the study participants are measured repeatedly
(at least twice) throughout the duration of the study, thereby permitting the
direct assessment of changes in the response variable over time [16]. To
illustrate, participants in a cross-sectional study likely vary in age;
however, this type of design cannot be used to study the effect of aging
because the effect of aging is potentially confounded with cohort effects. It
is important to note here that although we are sampling the evolving network
at multiple time points, we are only asking participant information that can
maintain privacy.

One example of a longitudinal network study is the work by Birkett et al.
[11]. The authors studied the impact of network-level mechanisms and STI
infections on the spread of HIV and found that network-level mechanisms and
STI infections play a significant role in the spread of HIV and in racial
disparities among (YMSM). Their work shows HIV prevention efforts should
target YMSM across race, and interventions focusing on YMSM partnerships with
older MSM might be highly effective. In general, one would expect observing a
network multiple times to provide more information, and therefore improve
accuracy of inference, compared to observing the network just once. In
addition, one may address questions that can only be interrogated in a
longitudinal study. When growing a network in a simulation, we can track every
iteration of the dynamic network and have arbitrarily many observations at our
disposal. In an actual study, one is of course constrained by resources and
logistics. If the data are obtained from self-administered or staff-
administered surveys, too frequent reporting may lead to participant burden
and reduce his or her willingness to continue participation, whereas too
infrequent reporting may lead to recall bias and participants may be lost to
follow up. For example, a person may not remember each individual whom they
dated over a five-year period and may not be able to reliably recall the
timing of the relationships. Collecting data at different time points that are
optimally spaced helps alleviate recall bias while still maintaining an avenue
for accurate inference.

Our goal in this paper is to implement a discrete-time version of the
mechanistic network model introduced by Hansson et al. and use the model to
identify optimal spacing between two data collection points (waves) in a
longitudinal network study such that we can achieve the dual goal of accurate
inference (learning model parameters as precisely as possible) while
minimizing participant burden (using network features that in practice could
be elicited from participants with a minimal number of survey questions).
These results have implications for study design for HIV and other sexually
transmitted diseases, and more broadly they can inform other research
questions involving (longitudinal) network data.

This paper is structured as follows. We discuss the discrete-time mechanistic
network model in Section 2.1 and explain our ABC-based approach to approximate
parameter inference in Section 2.2. We show our results in Section 3 and
conclude with a discussion in Section 4.

##  2 Methods

###  2.1 Mechanistic network model

As noted in the Introduction, there are several mechanistic models for MSM
contact formation in specific populations. We focus on the mechanistic model
introduced to study MSM contact networks in Stockholm, Sweden [13]. The model
incorporates specific behaviors that guide the formation and dissolution of
sexual contacts as well as migration of individuals in and out of the
population. While the original model was formulated in continuous time, we
consider a discrete time version of the model. This means that rates in the
original formulation correspond to probabilities in ours. We note that as the
number of the potential discrete time events tends to infinity and the event
probabilities tend to zero, our formulation of the model converges to the
original. Throughout this paper, each discrete model time step iteration is
taken to correspond to one calendar month, and all events are recorded at the
end of each iteration. While a constant number of individuals enter the
population at each iteration, each individual leaves the population with a
fixed probability at each iteration. The size of the network therefore
fluctuates around $n$ nodes, where $n$ is the initial number of nodes in the
network.

The model incorporates two types of partnerships: steady and casual. Casual
relationships are defined to only last one iteration at onset while steady
relationships are defined to have the potential to last longer. An individual
can have at most one steady partner at any given time. The probability that a
single person enters a steady relationship at a given iteration is $\rho
P_{0}$, where $P_{0}$ is the proportion of single individuals in the present
iteration. In the original model, where $\rho$ is a rate of steady partnership
formation, $P_{0}$ fluctuates around an equilibrium; in our version, we fix
this parameter and absorb it into $\rho$ for simplicity and to improve
identification of model parameters. Our modified probability of a single
person entering a steady relationship at a given iteration is therefore
$\rho$. While the number of people willing to form relationships varies from
iteration to iteration, the probability a single person joining a relationship
stays the same. In the Hansson paper, the differential equation formulation of
the model explicitly considers the fluctuation of the likelihood of new
relationships while we do it implicitly as the number of singles changes. The
probability of said steady relationship dissolving at each iteration is
$\sigma$.

In addition to a steady relationship, an individual may also have one casual
partnership at each iteration. These casual relationships may occur alongside
steady partnerships or during times when the person is single. A single
individual enters a casual relationship with probability $\omega_{0}$ while an
individual who is currently in a steady relationship forms a casual
relationship with probability $\omega_{1}$. For any partnership to form, both
individuals must be willing to join that relationship. In the scenario where
an odd number of individuals would like to form a relationship, one of them
(chosen at random) is left out. Each person migrates from the population with
probability $\mu$, and individuals enter into the population at constant rate
$n\mu$. The migration of an individual and the formation and dissolution of a
sexual contacts are all determined by the outcome of independent Bernoulli
trials. In the original continuous time formulation of the model, duration of
steady relationships and the time spent in the population both follow
exponential distributions. In contrast, for our discrete time formulation both
are geometrically distributed. Starting from an empty graph with $n$ nodes, we
first run the model until we are confident that it has converged to the target
distribution. We set our migration probability to 0 to ensure we are sampling
individuals longitudinally and to maintain a closed cohort design. We note
that the ’constant’ number of nodes being added is largely dependent on only
$\mu$ and $n$ and easily recoverable. The model is described in Algorithm 1
and a few graph realizations from the model are illustrated in Figure 1. We
chose 1000 iterations to ensure we are past the burn-in [13].

Figure 1: Network visualizations containing cumulative (from iteration 1)
steady (red dashed) and casual (blue solid) edges for iterations 1 (left), 6
(middle), and 12 (right). We used the following parameter values: $\mu$ = 0,
$\rho$ = 0.3, $\sigma$ = 0.1, $w_{1}$ = 0.2, $w_{0}$ = 0.4. Algorithm 1
Hansson MSM model [13]

1:Inputs:

2:$n$ := number of nodes

3:G := (V, E) the graph has no edges

4:$\rho$ := scaling parameter for partnership formation probability

5:$\sigma$ := separation probability for steady relationships

6:$\omega_{0}$ := the probability someone who is single enters into a casual
relationship

7:$\omega_{1}$ := the probability someone who is in a partnership enters into
a casual relationship

8:$iterations$ := the number of iterations to run algorithm

9:Algorithm:

10:for i in 1:iterations do

11: Dissolve all casual relationships generated from the previous iteration if
applicable

12: Identify all current steady relationships

13: if there are a positive number of steady relationships then

14: Dissolve each with probability $\sigma$

15: end if

16: Identify all nodes with degree 0

17: Set willingness to form a steady relationship with probability $\rho$

18: Randomly match the maximum number of even nodes that are willing to form a
steady relationship into edges

19: Identify single nodes := nodes with degree 0

20: Set willingness to form a casual relationship with probability
$\omega_{0}$

21: Identify steady nodes := nodes with degree 1

22: Set willingness to form a casual relationship with probability
$\omega_{1}$

23: Randomly match the maximum number of even nodes among the single and
steady nodes together that are willing to form a casual relationship

24: Record Network

25:end for

###  2.2 Inference of model parameters

In Bayesian inference, complete knowledge of the model parameters, given the
observed data, is contained in the posterior distribution. Typically, in
mechanistic models, the complexity of the model means that the likelihood and
corresponding posterior distribution is not available in closed form. In
mechanistic models one can nevertheless forward simulate data from the model
given parameters, and these parameter values may be obtained from a prior
distribution. ABC is an inference framework that has been developed to deal
with models that have intractable likelihoods. There are several ABC methods
to generate samples from an approximate posterior distribution. For clarity of
our objective, we use the simple accept/reject algorithm. In accept/reject, we
propose parameter values from a prior distribution to generate data and retain
parameter values that produce data that resembles the observed data [17]. If
we only kept parameter values that reproduced the observed data exactly, this
approach would recover the exact posterior for discrete data. This approach is
however impossible for continuous data because the probability of sampling a
continuous value exactly is 0. To ensure that our prior distribution’s support
is realistic to MSM relationship characteristics, we utilize a uniform
distribution on the duration of average time spent for people to be open to
joining a relationship $\frac{1}{\rho}$ [1 month, 50 months], average time a
steady relationship lasts $\frac{1}{\sigma}$ [1 month, 90 months], average
time for a single individual to partake in a casual relationship
$\frac{1}{\omega_{0}}$ [1 month, 40 months], and average time for an
individual in a relationship to partake in a casual relationship
$\frac{1}{\omega_{1}}$ [1 month, 61 months]. Figure 2 shows the prior
distributions on these inverse parameters and the corresponding implied prior
distributions on the parameters themselves. We recognize a variety of
definitions for steady and casual relationships in MSM contact networks, as
well as a variety of estimates for the support of each duration [18, 19, 20,
21, 22, 23, 24, 25]. We chose our support to be consistent with the data the
model was originally trained on [13], and calculate the reciprocal of each
parameter sampled from the prior as an input to our model.

Figure 2: Prior distributions on the inverse parameters and the corresponding
implied prior distributions on parameters themselves. The top row shows the
distributions of the inverse parameters, which can be interpreted as
distributions of the average values of geometric distributions. The bottom row
shows the distributions of the parameter values themselves for our discrete
time mechanistic network model for the following parameters: $\rho$, $\sigma$,
$\omega_{1}$, $\omega_{0}$.

There are at least three major considerations in the ABC accept/reject
framework: summary statistics, distance measure, and similarity threshold
[26]. Given the mechanistic network model of interest, we manually chose a set
of summary statistics needed for inference, which renders the network space
more manageable [26]. The choice of network summary statistics was guided by
the principle that it should be possible to obtain this information from study
participants using a questionnaire and they should be informative of the model
parameters. At a minimum, one needs at least as many summary statistics as
there are parameters to be inferred [26]. Although the model has five
parameters (six if one counts $n$), as previously mentioned, we opted to fix
one of them, the migration probability $\mu=0$. This leaves us with four
parameters to recover: probability of a single person entering a steady
relationship $\rho$, probability of dissolving a steady relationship $\sigma$,
probability of a single individual to enter a casual relationship
$\omega_{0}$, and probability of an individual in steady relationship to enter
a casual relationship $\omega_{1}$. More information on parameters can be
found in Table 1. We chose the four summaries, denoted $s_{1}$ through
$s_{4}$, as listed in Table 2. We chose summaries that could be elicited by
asking participants to consider their sexual history in the past year only.
Longer histories could potentially be more informative, but longer look-back
periods would likely increase recall bias.

Since the mechanistic network model is outside the exponential family, we have
no guarantee of sufficiency, i.e., that our summary statistics fully summarize
our network. However, we still require the summary statistics to be
informative of the model parameters. An informal way to assess the extent of
informativeness is to investigate plots of network summary statistics against
model parameters. We denote these relationships as $s_{i}(\theta)$ where
$i\in\\{1,2,3,4\\}$ and $\theta\in\\{\rho,\sigma,w_{0},w_{1}\\}$. We refer to
these relationships as _mapping functions_ , and we estimate them using
simulations where $s_{i,k}(\theta)$ represents the value of summary statistic
$s_{i}$ with respect to generative parameter $\theta$ in simulation run $k$.
The value of $s_{i}(\theta)$ is given as the median value of $s_{i,k}(\theta)$
taken across all simulations $k$.

We measured the distance between a simulated network and the observed network
by calculating the Euclidean distance within the normalized summary statistic
space. The normalized summary statistic value is obtained by first subtracting
the mean of the summary statistic from each value and then dividing each value
by the standard deviation of the summary statistic. We populated an ABC
reference table for each lag by generating 10,000 graphs by sampling
parameters from their joint priors and varying the lag between the two network
observations between zero and 150 iterations. Next, taking a sample per lag
from our joint prior density and its corresponding graph as our ground truth,
we simulated samples from the corresponding approximate posterior
distribution. We retained the parameters associated with the 100 (top 1%)
smallest distances in the normalized summary statistic space between the
observed and generated graphs.

Finally, we performed a regression adjustment on samples from the approximate
posterior distribution [27, 28]. The goal of the regression adjustment is to
improve our ABC posterior’s convergence to our target posterior. The basis of
the method is that we can obtain an estimate of our expected parameter values
given the summaries using linear regression in the localized neighborhood
around our observed data that we get from the approximate posterior. Then, we
can use this relationship to adjust our approximate posterior distribution
[27]. We normalized the parameters in our reference table, and utilized the
root mean squared error (RMSE) of the approximated posteriors for a fixed set
of 500 ground truth parameters to measure accuracy of inference. Then, we
averaged over all 500 parameter sets for an estimate of the RMSE, for a given
lag, over our prior space [29]. For clarity, consider $\theta_{i}$ as the
$ith$ ground truth parameter and $\hat{\theta}_{i,k}$ as the $kth$ sample from
the approximate posterior estimating $\theta_{i}$. Our estimate of RMSE is
then

$\hat{RMSE}=\frac{1}{500}\Sigma_{i=1}^{500}\sqrt{\frac{1}{100}\Sigma_{k=1}^{100}(\theta_{i}-\hat{\theta}_{i,k})^{2}}.$
(1)

Finally, we fitted a locally weighted regression (loess) with a 95% confidence
interval to the data. We note that while the values of the parameters in the
reference table are normalized, the resulting approximate posterior
distributions are not. Individually normalizing the reference table parameters
is useful because it places all parameters on the same scale when calculating
the RMSE. However, the posteriors are displayed on the original scale for ease
of interpretation.

Parameter |  Hansson et. al [13] |  Current Paper
---|---|---
n |  Average population size |  Population size
$\mu$ |  Rate of leaving the sexually active population |  Probability of leaving the sexually active population
$\rho$ |  Partnership desire scale rate |  Probability of a single person entering a steady relationship
$\sigma$ |  Separation rate |  Probability of dissolving a steady relationship
$\omega_{0}$ |  Rate at which an individual who is single tries to have casual sex |  Probability of a single individual to enter a casual relationship
$\omega_{1}$ |  Rate at which an individual who is in a relationship tries to have casual sex |  Probability of an individual in a steady relationship to enter a casual relationship
Table 1: Parameter And Their Interpretation Name |  Description
---|---
$s_{1}$ |  The proportion of single individuals
$s_{2}$ |  The average length of steady relationships that start and end in the course of the study
$s_{3}$ |  The proportion of individuals in steady relationships who are also in casual relationships
$s_{4}$ |  The proportion of steady relationships among all relationships
Table 2: Summary Statistic Descriptions

##  3 Results

We evaluated the mapping functions on a grid along the unit interval by
generating 100 graphs per parameter value and using box plots to summarize the
results. We investigated mapping functions in two different scenarios. First,
we varied each parameter in turn while keeping all others fixed at the values
reported in [13], i.e., we fixed $\rho$ = 0.3, $\sigma$ = 0.1, $\omega_{0}$ =
0.4, and $\omega_{1}$ = 0.2, and we also set $\mu$ = 0 to keep the cohort
closed. Second, we sampled each free parameter from its respective prior
distribution. These plots were used to ensure that the chosen summary
statistics are informative about the model parameters as can be seen in
Figures 3 and 4.

While more summaries could be included, that would increase the computational
burden and likely would not significantly increase accuracy. We considered
several extra summaries during discovery. Since we would like to obtain the
data from questionnaires, one also needs to consider participant burden: all
else equal, we would like to ask as few questions as needed to address the
scientific question at hand. It is also worth emphasizing that each of the
listed summaries can be obtained using privacy preserving questions only in
data collection, i.e., participants do not need to disclose their identity nor
the identity of their steady or casual partners. This arguably improves the
quality of the collected data as respondents would be expected to be more
likely to report their behavior accurately. We note that while a regression
adjustment generally improves the results, it can at times generate
functionally impossible values, such as negative probabilities, or worsen our
inference when the summaries do not accurately represent the network. In the
rare occasion the regression adjustment proposes a negative number, we opt to
take a conservative approach and set the value at 0. In this study, we did not
see any adjusted proposal probabilities above 1.

We visualized the regression adjusted approximate posteriors when looking at
the graph once or twice with a lag of 50 iterations in Figure 5. Furthermore,
as expected, and as shown in Figure 6, observing a network twice results in a
smaller average error compared to observing a network only once. The
improvement is largely driven by our ability to recover $\sigma$ and $\rho$
parameters as shown in Figure 7. We also see the average error steadily
decreases with the lag between the two network observations until about 40 to
50 iterations. This lag between the two network observations (data collection
waves) is optimal in the sense that extending the gap further does not greatly
increase accuracy of inference but does lengthen the duration of the study. In
a closed cohort study, all else equal, the longer the duration of the study,
the greater the expected attrition of study participants. Attrition of study
participants in a setting like ours would lead to incomplete ascertainment of
network structure and therefore introduce an additional error to network
summary statistics. We also see that implementing a regression adjustment does
reduce our average error by nearly an additional 2.6%, while maintaining the
overall trend and optimal lag. Finally, we note our overall ability to discern
parameters from our joint prior distribution when collecting data twice after
a regression adjustment with an optimal drop of roughly 62% from our average
prior error and 18% when only collecting data once.

Figure 3: Pairwise relationships between the model parameters (horizontal
axes) and the summary statistics (vertical axes) used in our ABC inference
scheme. Free parameters are fixed at $\mu$ = 0, $\rho$ = 0.3, $\sigma$ = 0.1,
$\omega_{0}$ = 0.4, $\omega_{1}$ = 0.2. The lag between two consecutive
network observations is fixed at 15 iterations. Each box plot consists of 100
samples. Figure 4: Pairwise relationships between the model parameters
(horizontal axes) and the summary statistics (vertical axes) used in our ABC
inference scheme. Free parameters are sampled from the prior distributions.
The lag between two consecutive network observations is fixed at 15
iterations. Each box plot consists of 100 samples. Figure 5: Approximate
marginal posterior distributions of model parameters obtained by retaining the
top 1% of proposed prior samples in our ABC accept/reject inference scheme.
Different rows correspond to comparing the prior (top), observing the graph
once (middle), and observing the graph twice with a lag of 50 iterations
(bottom). All posteriors include a regression adjustment. The blue solid lines
represent the 95% credible intervals and the red dotted lines represent the
true parameter values. Figure 6: Estimated average RMSE, where the average is
taken across multiple network realizations, as a function of the lag between
the two network observations. We also include a loess curve with a 95%
confidence interval (shaded areas). The average prior average RMSE is $2.22$
(not shown), whereas the corresponding regression adjusted error for a network
observed only once is $1.03$ that for a network observed twice with a lag of
50 iterations is $0.84$. Figure 7: Estimated regression adjusted average RMSE
for the total error (top curve) and separately for the four parameters
considered in our study (bottom four curves). These results show that when
observing a network twice, the reduction in total RMSE is mainly due to the
reduction of RMSE for $\rho$ and $\sigma$.

##  4 Discussion

In this paper, we investigated the accuracy of an approximate inference scheme
applied to an evolving mechanistic network model in a setting where the
network, representing sexual contacts among people in a closed population, is
observed at two different time points. As expected, observing the network
twice improves the accuracy of inference, but this reduction in inferential
error depends on the time lag between the two observations. Given that
collection of real-world sexual network data is expensive and logistically
challenging, it pays off to optimize the gap between the two time points to
maximize accuracy of inference. If the two network observations are too close
in time, there may have been only minimal changes in the network structure,
and therefore the second observation adds little information. However, if the
two network observations are too far apart in time, the study may be
logistically difficult to carry out in practice and the population is likely
to experience significant churn.

There are a total of six parameters in the model, but we fixed two of them to
focus on a closed, fixed-sized cohort. When considering the contribution of
the remaining four parameters to inferential error, we observed that the
$\sigma$ (probability of dissolving a steady relationship) and $\rho$
(probability of a single person entering a steady relationship) parameters
benefited the most from the lag between the two network observations. This
finding is intuitive as these two parameters influence multiple relationship
iterations. However, $\omega_{0}$ (probability of a single individual to enter
a casual relationship) and $\omega_{1}$ (probability of an individual in
steady relationship to enter a casual relationship) both correspond to one-
time events and do not benefit as much from a lag. In particular, $\omega_{1}$
is relatively accurate at all lags while $\omega_{0}$ would likely see more
relative improvement through the consideration of another summary statistic.

The set of summary statistics that may be considered in inference depends on
the information obtained from subjects through study questionnaires. The
informativeness of questions themselves depends on the mechanisms that drive
contact formation in the study population. Depending on the mechanisms, it is
possible that any set of individual-level questions (giving rise to so-called
egocentric samples of the network) may be inadequate for network inference and
instead one may need information about the full network structure. While this
type of network-level information could be obtained using a sociocentric
design, it is very challenging, and we are aware of only one study that has
implemented this in practice. The Likoma Network Study was based on a
sociocentric survey of sexual partnerships aimed to investigate the
population-level structure of sexual networks connecting the young adult
population of several villages on Likoma Island, Malawi [30]. We stress that
this notable study is cross-sectional and therefore corresponds to a one-time
observation of the network (even if the data collection in this study occurred
in two stages for logistical reasons). Obtaining two observations of the
network would be logistically nearly impossible, and doing so in larger
populations is not feasible.

Our results highlight the importance of using simulation to investigate the
hypothesized generative mechanisms of network formation to inform future study
designs, here specifically 1) what questions to ask so that maximally
informative network summary statistics may be constructed and 2) how to space
the two (or possibly more) data collection waves. For example, in our setting,
introducing extensive migration in the population leads to a shorter optimal
lag between the two network observations. Our approach is compatible with the
recommended paradigm of using simulations for designing and interpreting
intervention trials in infectious diseases, particularly with regard to
emerging infectious diseases [31]. One of the main goal of such simulations is
to more accurately reflect the dynamics of the transmission process. For
sexually transmitted diseases, learning about the mechanisms of network
formation is an important step in that direction.

In this paper, we have used basic ABC and basic regression adjustment
techniques because our goal here is to see whether the ABC approach is
effective in its simplest and most interpretable form. More refined variants
of these methods, which can substantially improve computational performance,
can be studied later on. Finally, at the time of writing, we came across
related work on how design choices for egocentric network studies impact
statistical estimation and inference for ERGMs [32]. This investigation is
relevant for ours, although our focus is specifically on the multiple
observation of the evolving network. For a suitably chosen ERGM, i.e., an ERGM
with reasonably simple dependence assumptions, it is possible to attain
sufficient summary statistics from egocentric network samples. This allows for
exact statistical inference, but at the cost of making distributional
assumptions that may not hold. For that reason, it is valuable for
investigators to have various methods at their disposal so that they may
choose the tool that best fits the scientific problem at hand.

##  5 Declarations

###  5.1 Availability of data and materials

Data was simulated using a mechanistic model introduced to study MSM contact
networks in Stockholm, Sweden [13]. Our code is accessible at:
https://github.com/onnela-lab/longitudinal-inference

###  5.2 Competing Interests

The authors have no competing interest to report.

###  5.3 Funding

NIH Award #R01AI138901.

###  5.4 Authors’ contributions

All authors conceived the study as well as drafted and revised the manuscript;
OS implemented the method in code and carried out data analyses. TH and JP
supervised.

###  5.5 Acknowledgments

We would like to acknowledge John Quackenbush and Marcello Pagano for their
thoughtful insights on summary statistic exploration.

## References

  *   * Macal et al. [2004] Macal, C., Sallach, D., North, M.: Emergent structures from trust relationships in supply chains. In: Proc. Agent 2004: Conf. on Social Dynamics, pp. 7–9 (2004)
  * Scholtens and Gentleman [2005] Scholtens, D., Gentleman, R.: Making sense of high-throughput protein-protein interaction data. Statistical Applications in Genetics and Molecular Biology 3(1) (2005)
  * Le et al. [2022] Le, T.-M., Raynal, L., Talbot, O., Hambridge, H., Drovandi, C., Mira, A., Mengersen, K., Onnela, J.-P.: Framework for assessing and easing global COVID-19 travel restrictions. Scientific Reports 12(1), 1–13 (2022)
  * Adamic and Huberman [2000] Adamic, L.A., Huberman, B.A.: Power-law distribution of the world wide web. Science 287(5461), 2115–2115 (2000)
  * Robins et al. [2007] Robins, G., Pattison, P., Kalish, Y., Lusher, D.: An introduction to exponential random graph (p*) models for social networks. Social Networks 29(2), 173–191 (2007)
  * Goyal and Onnela [2020] Goyal, R., Onnela, J.: Framework for converting mechanistic network models to probabilistic models. arXiv 2001.08521 (2020)
  * Albert and Barabási [2002] Albert, R., Barabási, A.-L.: Statistical mechanics of complex networks. Reviews of Modern Physics 74(1), 47 (2002)
  * Wertheim et al. [2011] Wertheim, J.O., Kosakovsky Pond, S.L., Little, S.J., De Gruttola, V.: Using HIV transmission networks to investigate community effects in HIV prevention trials. PloS One 6(11), 27775 (2011)
  * Aroke et al. [2022] Aroke, H., Katenka, N., Kogut, S., Buchanan, A.: Network-based analysis of prescription opioids dispensing using exponential random graph models (ERGMs). In: Complex Networks & Their Applications X: Volume 2, Proceedings of the Tenth International Conference on Complex Networks and Their Applications COMPLEX NETWORKS 2021 10, pp. 716–730 (2022). Springer
  * Rolls et al. [2013] Rolls, D.A., Wang, P., Jenkinson, R., Pattison, P.E., Robins, G.L., Sacks-Davis, R., Daraganova, G., Hellard, M., McBryde, E.: Modelling a disease-relevant contact network of people who inject drugs. Social Networks 35(4), 699–710 (2013)
  * Birkett et al. [2015] Birkett, M., Armbruster, B., Mustanski, B., et al.: A data-driven simulation of HIV spread among young men who have sex with men: the role of age and race mixing, and STIs. Journal of Acquired Immune Deficiency Syndromes 70(2), 186 (2015)
  * Mei et al. [2010] Mei, S., Sloot, P.M., Quax, R., Zhu, Y., Wang, W.: Complex agent networks explaining the HIV epidemic among homosexual men in Amsterdam. Mathematics and Computers in Simulation 80(5), 1018–1030 (2010)
  * Hansson et al. [2019] Hansson, D., Leung, K.Y., Britton, T., Strömdahl, S.: A dynamic network model to disentangle the roles of steady and casual partners for HIV transmission among MSM. Epidemics 27, 66–76 (2019)
  * Padeniya [2021] Padeniya, S.M.T.N.: Mathematical modelling to explore the role of the female-sex-worker-client interaction for gonorrhoea transmission and prevention among australian heterosexuals. PhD thesis, UNSW Sydney (2021)
  * Vajdi et al. [2020] Vajdi, A., Juher, D., Saldaña, J., Scoglio, C.: A multilayer temporal network model for STD spreading accounting for permanent and casual partners. Scientific Reports 10(1), 1–12 (2020)
  * Fitzmaurice et al. [2012] Fitzmaurice, G.M., Laird, N.M., Ware, J.H.: Applied longitudinal analysis (2012)
  * Csilléry et al. [2010] Csilléry, K., Blum, M.G., Gaggiotti, O.E., François, O.: Approximate bayesian computation (ABC) in practice. Trends in Ecology & Evolution 25(7), 410–418 (2010)
  * Malone et al. [2018] Malone, J., Syvertsen, J.L., Johnson, B.E., Mimiaga, M.J., Mayer, K.H., Bazzi, A.R.: Negotiating sexual safety in the era of biomedical HIV prevention: relationship dynamics among male couples using pre-exposure prophylaxis. Culture, Health & Sexuality 20(6), 658–672 (2018)
  * Down et al. [2017] Down, I., Ellard, J., Bavinton, B.R., Brown, G., Prestage, G.: In Australia, most HIV infections among gay and bisexual men are attributable to sex with ‘new’partners. AIDS and Behavior 21(8), 2543–2550 (2017)
  * de Vroome et al. [2000] Vroome, E.M., Stroebe, W., Sandfort, T.G., WIT, J.B., Griensven, G.J.: Safer sex in social context: Individualistic and relational determinants of AIDS-preventive behavior among gay men 1. Journal of Applied Social Psychology 30(11), 2322–2340 (2000)
  * Wall et al. [2013] Wall, K.M., Stephenson, R., Sullivan, P.S.: Frequency of sexual activity with most recent male partner among young, internet-using men who have sex with men in the United States. Journal of Homosexuality 60(10), 1520–1538 (2013)
  * Davidovich [2006] Davidovich, E.: Liaisons dangereuses: HIV risk behavior and prevention in steady gay relationships (2006)
  * Weiss et al. [2020] Weiss, K.M., Goodreau, S.M., Morris, M., Prasad, P., Ramaraju, R., Sanchez, T., Jenness, S.M.: Egocentric sexual networks of men who have sex with men in the United States: Results from the ARTnet study. Epidemics 30, 100386 (2020)
  * Bavinton et al. [2016] Bavinton, B.R., Duncan, D., Grierson, J., Zablotska, I.B., Down, I.A., Grulich, A.E., Prestage, G.P.: The meaning of ‘regular partner’in HIV research among gay and bisexual men: implications of an australian cross-sectional survey. AIDS and Behavior 20(8), 1777–1784 (2016)
  * Myers et al. [1999] Myers, T., Allman, D., Calzavara, L., Morrison, K., Marchand, R., Major, C.: Gay and bisexual men’s sexual partnerships and variations in risk behaviour (1999)
  * Sisson et al. [2018] Sisson, S.A., Fan, Y., Beaumont, M.: Handbook of approximate bayesian computation (2018)
  * Beaumont [2019] Beaumont, M.A.: Approximate bayesian computation. Annual Review of Statistics and its Application 6, 379–403 (2019)
  * Beaumont et al. [2002] Beaumont, M.A., Zhang, W., Balding, D.J.: Approximate bayesian computation in population genetics. Genetics 162(4), 2025–2035 (2002)
  * Fearnhead and Prangle [2012] Fearnhead, P., Prangle, D.: Constructing summary statistics for approximate bayesian computation: semi-automatic approximate bayesian computation. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74(3), 419–474 (2012)
  * Helleringer and Kohler [2007] Helleringer, S., Kohler, H.-P.: Sexual network structure and the spread of HIV in Africa: evidence from Likoma Island, Malawi. Aids 21(17), 2323–2332 (2007)
  * Halloran et al. [2017] Halloran, M.E., Auranen, K., Baird, S., Basta, N.E., Bellan, S.E., Brookmeyer, R., Cooper, B.S., DeGruttola, V., Hughes, J.P., Lessler, J., et al.: Simulations for designing and interpreting intervention trials in infectious diseases. BMC Medicine 15(1), 1–8 (2017)
  * Krivitsky et al. [2022] Krivitsky, P.N., Morris, M., Bojanowski, M.: Impact of survey design on estimation of exponential-family random graph models from egocentrically-sampled data. Social Networks 69, 22–34 (2022)


# Congeniality Bounds on Quark Masses from Nucleosynthesis

M. Hossain Ali<EMAIL_ADDRESS>Department of Applied Mathematics,
Rajshahi University, Rajshahi 6205, Bangladesh  M. Jakir Hossain
<EMAIL_ADDRESS>Department of Mathematics, Rajshahi University, Rajshahi
6205, Bangladesh  Abdullah Shams Bin Tariq<EMAIL_ADDRESS>Department of
Physics, Rajshahi University, Rajshahi 6205, Bangladesh

###### Abstract

The work of Jaffe, Jenkins and Kimchi [Phys. Rev. D79, 065014 (2009)] is
revisited to see if indeed the region of congeniality found in their analysis
survives further restrictions from nucleosynthesis. It is observed that much
of their congenial region disappears when imposing conditions required to
produce the correct and required abundances of the primordial elements as well
as ensure that stars can continue to burn hydrogen nuclei to form helium as
the first step in forming heavier elements in stellar nucleosynthesis. The
remaining region is a very narrow slit reduced in width from around 29 MeV
found by Jaffe et al. to only about 2.2 MeV in the difference of the
nucleon/quark masses. Further bounds on $\delta m_{q}/m_{q}$ seem to reduce
even this narrow slit to the physical point itself.

###### pacs:

14.65.Bt, 26.20.Cd, 98.80.Ft

A reasonably contemporary approach is to study, even without going into
anthropic arguments, the nature of alternative universes as one changes the
values of physical parameters. In the parameter space, one then looks for
regions that could be similar to our universe and may possibly be congenial to
the creation and sustenance of intelligent life Weinberg (1989); Tegmark
(1997, 1998); Tegmark and Rees (1998); Hogan (2000); Tegmark  _et al._ (2006);
Jaffe  _et al._ (2009); Donoghue  _et al._ (2010). Bounds thus obtained may be
referred to as congeniality bounds.

In a recent work, Jaffe, Jenkins and Kimchi Jaffe  _et al._ (2009) studied how
sensitive our universe would be to variations of quark masses. For this they
chose to study the variations of masses of the three lightest quarks $u$, $d$
and $s$, under the constraint that the sum of these masses, $m_{T}$ remained
fixed. They also studied variations of $m_{T}$.

Their basic idea was to find the two lightest baryons for any quark mass
combination and consider them to play the roles of the proton and neutron in
forming nuclei. In this process they also considered $\Lambda_{\mathrm{QCD}}$
to be an adjusted free parameter that they tuned to keep the average nucleon
mass at 940 MeV. They then studied the variation of nuclear stability and, in
light of this, tried to obtain the regions of the parameter space where
nuclear chemistry in a somewhat familiar form could be sustained.

The starting point is that the three light quark masses would be changed
keeping their sum $m_{T}$ fixed. This parameter space can be neatly shown in
the form of an equilateral triangle [Fig. 1] where the distances of a point
from the base and right and left sides are, respectively, the masses of the
up, down and strange quarks.

Figure 1: The model space of light quark masses for a fixed $m_{T}$ shown in
the form of a triangle where the distance from the three sides give the three
masses. Figure reproduced from Jaffe  _et al._ (2009).

In this manner they identified congenial regions in a triangle parametrized in
terms of $x_{3}$ and $x_{8}$ [Fig. 2] defined as

$\displaystyle x_{3}$ $\displaystyle=$
$\displaystyle\frac{2m_{3}}{\sqrt{3}}\frac{100}{m_{T}^{\oplus}}=\frac{100\left(m_{u}-m_{d}\right)}{\sqrt{3}m_{T}^{\oplus}}$
(1) $\displaystyle x_{8}$ $\displaystyle=$
$\displaystyle\frac{2m_{8}}{\sqrt{3}}\frac{100}{m_{T}^{\oplus}}=\frac{100\left(m_{u}+m_{d}-2m_{s}\right)}{\sqrt{3}m_{T}^{\oplus}}$
(2)

Figure 2: The model space of light quark masses parametrized in terms of
$x_{3}$ and $x_{8}$ reproduced from Jaffe  _et al._ (2009). The point labeled
‘us’ points to the physical value in our present universe and therefore has
coordinates $(x_{3}^{\oplus},x_{8}^{\oplus})$.

Here $m_{T}^{\oplus}$ is the sum of the light ($u$, $d$ and $s$) quark masses
in our universe. It is obvious that $x_{3}$ basically gives the isospin
splitting while $x_{8}$ is related to the breaking of $SU(3)_{f}$ due to the
mass of the strange quark. Their results can be summarized in Fig. 3, where
the congenial regions are indicated in green 111If you are reading a black and
white print, then green appears as lightly shaded and red as deep shaded. This
triangle is for $m_{T}^{\oplus}$, i.e. with $m_{T}$ as it is in the present
universe. They have also studied variations in $m_{T}$, but our work is
limited to commenting on the case of $m_{T}^{\oplus}$, understanding that the
same arguments qualitatively extend to other values of $m_{T}$. This is
further justified in the discussions near the end of this report.

Figure 3: (color online) Figure reproduced from Jaffe  _et al._ (2009)
identifying congenial regions in the quark mass triangle with green bands. The
red and white regions are uncongenial and uncertain, respectively.

It is increasingly being understood that if there is complexity, fine-tuning
is inevitable Bradford . Even if one is not happy with anthropic arguments, we
simply cannot get away from fine-tuning. With this in mind, the first
impression that one has from Fig. 3, is that the congenial region seems to be
surprisingly large - allowing, around one order of magnitude variations in the
quark masses. Though this already involves intricate compensating adjustments
in $\Lambda_{\mathrm{QCD}}$ to keep the average ‘nucleon’ mass fixed.

However, one should appreciate the difficulty in setting up a new framework in
which a problem can be studied. From this perspective the authors of Jaffe
_et al._ (2009) should be commended for presenting, literally from scratch, a
setup for studying the congeniality bounds on quark masses. This setup can be
extended removing some of the constraints used in any further work. Indeed,
considering the significance of the work it was chosen first for a Viewpoint
article in Physics Perez (2009) and then went into a cover story in Scientific
American Jenkins and Perez (2010).

In our work, we remain within the provided setup, but extend the analysis to
bounds provided by nucleosynthesis. It should be noted here that whereas, on
the one hand, nuclear masses and stability expectedly vary comparatively
slowly with quark masses; on the other hand, the observed abundances of the
lightest nuclei hydrogen and helium provide much more stringent bounds on the
variation of nucleon masses. We report below how the congeniality triangle of
Fig. 3 is modified by the application of these constraints.

At the outset, the variation of the octet baryon masses as one traverses along
the borders of the triangle (Figs. 10 and 11 of Jaffe  _et al._ (2009)) were
reproduced to gain confidence in our code and our understanding of the
framework. The fitted parameter values of $c_{T}$, $c_{3}$ and $c_{8}$ from
Table III of Ref. Jaffe  _et al._ (2009) were used in the equation

$M_{B}=C_{0}+c_{T}x_{T}+c_{8}x_{8}+c_{3}x_{3}+\left\langle
B\left|H_{\mathrm{EM}}\right|B\right\rangle.$ (3)

This leads to

$\displaystyle M_{p}=C_{0}+3.68\>x_{T}+3.53\>x_{8}+1.24\>x_{3}+0.63$ (4)
$\displaystyle M_{n}=C_{0}+3.68\>x_{T}+3.53\>x_{8}-1.24\>x_{3}-0.13.$ (5)

The quantity occurring most in the analysis below being

$M_{n}-M_{p}=-2.48\>x_{3}-0.76.$ (6)

It may be noted here that using updated values of baryon masses from the
Particle Data Book changes the parameters very slightly and this is neglected
considering the qualitative nature of this work. After a clarification on the
adjustment of $C_{0}$ (corresponding to an adjustment of
$\Lambda_{\mathrm{QCD}}$) from the authors Jaffe  _et al._ it was possible to
reproduce the figures.

Then the issue of further bounds from nucleosynthesis were studied. It is well
known that the observed abundances of the primordial nuclei hydrogen
(protons), helium (alpha particles) etc. are sensitively tied to the masses of
the nucleons Barrow and Tipler (1986); Hogan (2000). The slight difference in
the masses of the proton and neutron are responsible for the survival of
protons with the observed abundance. The (un)congenial regions of the triangle
is explored further under these constraints.

There are three cases that arise here:

### Case I: $x_{3}>x_{3}^{\oplus}$

Let us concentrate on the region on the upper right of the triangle with
$x_{3}$ values greater than at the point labeled ‘us’ on the right hand side
of the triangle.

Of the two nucleons, the neutron is heavier by about 1.3 MeV. If it was just
0.8 MeV less, that would bring it below the electron capture threshold for
protons, i.e. it would become energetically favourable for protons to capture
electrons and become neutrons. All the protons would have been converted to
neutrons in the Big Bang. The Universe would be full of neutrons and nothing
else. We would not be here. In words of Barrow and Tipler Barrow and Tipler
(1986)

> Without electrostatic forces to support them, solid bodies would collapse
> rapidly into neutron stars or black holes. Thus, the coincidence that allows
> protons to partake in nuclear reactions in the early universe also prevents
> them decaying by weak interactions. It also, of course, prevents the 75% of
> the Universe which emerges from nucelosynthesis in the form of protons from
> simply decaying away into neutrons. If that were to happen no atoms would
> ever have formed and we would not be here to know it. [Ref. Barrow and
> Tipler (1986), p. 400]

The same issue is also discussed by Hogan Hogan (2000)

> The $u$-$d$ mass difference in particular attracts attention because the $d$
> is just heavier enough than $u$ to overcome the electromagnetic energy
> difference to make the proton ($uud$) lighter than the neutron ($udd$) and
> therefore stable. On the other hand, if it were a little heavier still, the
> deuteron would be unstable and it would be difficult to assemble any nuclei
> heavier than hydrogen.

Therefore, it is necessary to have

$M_{n}-M_{p}\geq 0.5\>\mathrm{MeV}.$ (7)

This reduces the congenial region on the upper right of the physical point to

$x_{3}\leq-0.51.$ (8)

### Case II: $x_{3}<x_{3}^{\oplus}$

Now let us move to the bottem-left side of the triangle (left of the $x_{8}$
axis) where $x_{3}$ values are smaller than that at the ‘us’-labeled point,
again concentrating on the right hand side of the triangle.

The key reaction by which hydrogen ‘burns’ in stars such as the sun involves
the reaction

$\displaystyle p+p$ $\displaystyle\rightarrow$ $\displaystyle
d+e^{+}+\nu+0.42\,\mathrm{MeV}$ (9) $\displaystyle e^{+}+e^{-}$
$\displaystyle\rightarrow$ $\displaystyle 1\>\mathrm{MeV}$ (10)

So the total amount of energy released in this reaction is 1.42 MeV.

If the neutron mass was 1.42 MeV (0.15%) more than it is, this reaction would
not happen at all. It would need energy to make it go, rather than producing
energy. Deuterons are a key step in burning hydrogen to helium. Without them,
hydrogen would not burn, and there would be no long-lived stars and no stellar
nucleosynthesis to produce the remaining elements.

Therefore, it is necessary that

$M_{n}-M_{p}\leq 2.72\>\mathrm{MeV}.$ (11)

This reduces the congenial region on the lower left of the physical point to

$x_{3}\geq-1.4.$ (12)

These two conditions, thus, significantly reduce the congenial corridor from

$-12.9\leq x_{3}\leq 4.1$ (13)

to

$-1.4\leq x_{3}\leq-0.5.$ (14)

It may be noted here that the width of this region is of the same order as the
uncertainty in $x_{3}^{\oplus}$ itself due to uncertainties in the light quark
masses given by $x_{3}^{\oplus}=-1.17\pm 0.43$.

In fact it was a pleasant surprise to realise, rather late into our work, that
Hogan Hogan (2000) reached essentially similar conclusions which were
expressed in terms of the up-down quark mass difference, $\delta m_{d-u}$ and
Section IV of his review Hogan (2000) is a recommended read for anybody
interested in this issue. The approximately 1.4 + 0.8 = 2.2 MeV window of
variation that we find is in agreement with the allowed region in Fig. 1 of
the same Hogan (2000).

### Case III: Left half of the triangle

If we move to the left half of the triangle, we essentially replace the down
quark with a strange quark. We know that the $s$-quark is, in some ways, like
a heavy $d$-quark. In the left half of the triangle the $s$-quark is light and
the $d$-quark is heavy. As if they simply interchange positions. That is why
Jaffe _e_ t al. Jaffe  _et al._ (2009) seem to find a symmetric congenial
region in the left of the triangle. The discussions for Cases I and II narrow
it down, but do not remove it.

However, let us now turn towards the coupling between $u$-$d$ and $u$-$s$. The
$u$-$d$ coupling is much stronger, whereas the $u$-$s$ coupling is suppressed.
This is described by the well-known Cabibbo angle $\theta_{C}$. Where the
$u$-$d$ coupling carries a factor $cos\,\theta_{C}$ and the $u$-$s$ coupling
carries a factor of $sin\,\theta_{C}$, the Cabibbo angle being about 13
degrees.

This is like the present world with a much weaker weak interaction. This the
case where the weak decay rate of neutrons is not strong enough to produce the
primordial neutron-proton abundance ratio of 1:6. Without this we are left
without enough protons, i.e. without enough hydrogen, which is key to both
stellar burning and biological life itself. Therefore we are left with only a
narrow region on the right [Fig. 4].

Figure 4:  (color online) Fig. 3 adapted by the further restrictions imposed
leaving only a very narrow congenial slit in the bottom-right region of the
triangle.

The only remaining question is probably regarding the length of this narrow
region extending nearly up to the centre of the triangle. As one moves up this
narrow slit towards the centre, away from the physical point, the up-down
quarks become heavier keeping the down quark slightly heavier than the up.
Meanwhile the strange quark becomes lighter to keep $m_{T}$ fixed. The physics
considered here is probably not very sensitive to the strange quark mass. The
increase in the up-down masses is offset by the compensating adjustment in
$\Lambda_{\mathrm{QCD}}$ to keep the nucleon masses fixed. Therefore, the
length of this region could probably be an artifact of the simultaneous and
compensating tuning of quark masses and $\Lambda_{\mathrm{QCD}}$.

This indeed has been one of the conclusions in Jaffe  _et al._ (2009) as
summarized more elegantly in Jenkins and Perez (2010); as well as Bradford
where, reviewing the alternative universe landscapes studied by Aguirre
(2001); Harnik  _et al._ (2006); Adams (2008); Jaffe  _et al._ (2009); Jenkins
and Perez (2010) it has been observed that if one is prepared to adjust
another parameter in a compensating manner, it might be possible to find other
regions in the parameter space that are also congenial. However, that does not
remove the fine-tuning problem, as the alternative values are still finely
tuned and this is inevitable to produce complexity as observed in our present
universe. Here most of the alternatives are removed and the narrow region
remains as a result of the compensating adjustments of
$\Lambda_{\mathrm{QCD}}$.

Indeed along the narrow region the sum of the two lightest quarks vary with
the strange quark mass going in the opposite direction to keep $m_{T}$ fixed.
If the effect of this could be quantified, it would probably be possible to
restrict even the length of the narrow region [See the Addendum].

For example, as noted by Hogan Hogan (2000),

> … the sum of the (up and down) quark masses controls the pion mass, so
> changing them alters the range of the nuclear potential and significantly
> changes nuclear structure and energy levels. Even a small change radically
> alters the history of nuclear astrophysics, for example, by eliminating
> critical resonances of nucleosynthesis needed to produce abundant carbon.
> 222An interesting additional note is that, here Hogan cites Hoyle, F., D. N.
> F. Dunbar, W. A. Wenzel, and W. Whaling, 1953, Phys. Rev. 92, 1095. This has
> been cited several times in different papers, sometimes with Phys. Rev.
> Lett. as the source, and a few times with the title “A state in C12
> predicted from astrophysical evidence”. However, we failed to find any such
> article and would appreciate any information on this reference. The nearest
> match was D. N. F. Dunbar, R. E. Pixley, W. A. Wenzel, and W. Whaling, 1953,
> Phys. Rev. 92, 649, an article on the resonance in 12C often dubbed the
> ‘Hoyle resonance’.

Here it should be added that a more up-to-date view is that the strongest
effect on the scalar scattering lengths and deuteron binding energy seem to be
due to the “sigma-resonance” exchange (or correlated two-pion scalar-isoscalar
exchange) dependence on $m_{\pi}$ Hanhart  _et al._ (2008); Pelaez and Rios
(2010).

As mentioned at the outset, the analysis here has been limited to the case of
$m_{T}=m_{T}^{\oplus}$. It has been noted by Jaffe et al. Jaffe  _et al._
(2009) that the widths of the two major congenial bands on the bottom-left and
bottom-right of the triangle are independent of $m_{T}$. Therefore, naturally
the further exclusions for $x_{3}>x_{3}^{\oplus}$, $x_{3}<x_{3}^{\oplus}$
reducing the width of the band should also apply to other values of $m_{T}$.
The exclusion of the left half of the triangle should also extend to other
values of $m_{T}$. Therefore, in summary, it can be expected that for all
values of $m_{T}$, after applying constraints from nucleosynthesis, there will
only remain a similar very narrow congenial band at the bottom-right of the
triangle.

An additional comment is due here on the possibilities of universes with
deuterons, sigma-hydrogen, or delta-helium playing the roles of hydrogen as
listed in Jenkins and Perez (2010) as a summary of Jaffe  _et al._ (2009). The
point made here does not contradict that these could be stable lightest
elements. It is only pointed out that stability alone is not enough to produce
and sustain nuclear chemistry in a manner familiar to us. Correct primordial
abundances and conditions for sustained stellar burning provide constraints
that are much more difficult to satisfy. This probably calls for a closer
analysis of the other half in Jenkins and Perez (2010) related to possible
universes without any weak interaction of Harnik  _et al._ (2006) where there
indeed has been a detailed discussion of these issues pertaining to
nucleosynthesis. However, that would have to be another project; whereas, this
work is focused on Jaffe  _et al._ (2009).

In summary, it can be observed that, primordial nuclear abundances and
processes of stellar nucleosynthesis provide much more stringent constraints
on quark masses than nuclear stability. Using these constraints it is possible
to significantly reduce the congenial region in the space of light quark
masses.

## Addendum

Our attention has been drawn through referee comments to studies of the bounds
from nucleosynthesis Bedaque  _et al._ (2011); Berengut  _et al._ (2013), the
latter appearing after the initial submission of this paper, on $\delta
m_{q}/m_{q}$, where $m_{q}$ is the average of the light (up and down) quark
mass and $\delta m_{q}$ is the change in $m_{q}$ keeping $m_{u}/m_{d}$ fixed.
Coincidentally, along the length of the remaining narrow congenial region
$m_{u}/m_{d}$ is approximately constant. The latest value is $\left|\delta
m_{q}/m_{q}\right|<0.009$ Berengut  _et al._ (2013). There are other values in
literature, but they are generally of the same order. Let us try to do a crude
estimate of the effect of this constraint. For $m_{q}\approx 3.8$ MeV, $\delta
m_{q}\approx 0.035$ MeV.

From eqs. 8 and 9, we get

$M_{N}=\left(M_{n}+M_{p}\right)/2=C_{0}+3.68\>x_{T}+3.53\>x_{8}+2.5,$ (15)

which given that $x_{T}$ is kept fixed, leads to

$\delta M_{N}=\delta C_{0}+3.53\>\delta x_{8}.$ (16)

Now, $x_{8}$ as defined in eq. 2 can be re-expressed in terms of $m_{q}$ as

$x_{8}=\frac{200\left(m_{q}-m_{s}\right)}{\sqrt{3}m_{T}^{\oplus}}.$ (17)

If $m_{T}$ is kept fixed then $\delta m_{s}=-\delta m_{q}$, leading to

$\delta x_{8}=\frac{400\left(\delta
m_{q}\right)}{\sqrt{3}m_{T}^{\oplus}}=0.08,$ (18)

where we have used $\delta m_{q}\approx 0.035$ MeV and $m_{T}^{\oplus}\approx
100$ MeV. The remaining congenial region is too small to show on a figure of
this scale. In fact the region $x_{8}=x_{8}^{\oplus}\pm 0.08$ is too small to
show on a plot of this scale and is also very small compared to the
uncertainty in the value of $x_{8}^{\oplus}$ itself $x_{8}^{\oplus}=-59.5\pm
1.1$ due to the uncertainties in the determination of the light quark masses.
However, one should remember our estimate is rather crude without appropriate
consideration of the uncertainties in $m_{q}$. Taking these into account will
increase the region, but keep it within the same order as the uncertainty in
$x_{8}^{\oplus}$ itself. In short there is practically no congenial region
outside $\left(x_{3}^{\oplus},x_{8}^{\oplus}\right)$.

###### Acknowledgements.

The authors are grateful to Robert L. Jaffe, Alejandro Jenkins and Itamar
Kimchi for their kind replies to our queries and helpful suggestions. MHA and
ASBT would also like to acknowledge the support of the Abdus Salam
International Centre for Theoretical Physics (ICTP), Trieste, Italy through a
Regular Associateship and a Junior Associateship, respectively.

## References

  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * Tegmark (1997) M. Tegmark, Class.Quant.Grav. 14, L69 (1997), arXiv:gr-qc/9702052 [gr-qc] .
  * Tegmark (1998) M. Tegmark, Annals Phys. 270, 1 (1998), arXiv:gr-qc/9704009 [gr-qc] .
  * Tegmark and Rees (1998) M. Tegmark and M. J. Rees, Astrophys.J. 499, 526 (1998), arXiv:astro-ph/9709058 [astro-ph] .
  * Hogan (2000) C. J. Hogan, Rev. Mod. Phys. 72, 1149 (2000), arXiv:astro-ph/9909295 [astro-ph] .
  * Tegmark  _et al._ (2006) M. Tegmark, A. Aguirre, M. J. Rees, and F. Wilczek, Phys.Rev. D73, 023505 (2006), arXiv:astro-ph/0511774 [astro-ph] .
  * Jaffe  _et al._ (2009) R. L. Jaffe, A. Jenkins, and I. Kimchi, Phys. Rev. D79, 065014 (2009), arXiv:0809.1647 [hep-ph] .
  * Donoghue  _et al._ (2010) J. F. Donoghue, K. Dutta, A. Ross, and M. Tegmark, Phys.Rev. D81, 073003 (2010), arXiv:0903.1024 [hep-ph] .
  * Note (1) If you are reading a black and white print, then green appears as lightly shaded and red as deep shaded.
  * (10) R. Bradford, Int. J. Theor. Phys. 50, 1577.
  * Perez (2009) G. Perez, Physics 2, 21 (2009).
  * Jenkins and Perez (2010) A. Jenkins and G. Perez, Sci. Am. 2010, 302 (2010).
  * (13) R. L. Jaffe, A. Jenkins, and I. Kimchi, Private Communication.
  * Barrow and Tipler (1986) J. D. Barrow and F. J. Tipler,  _The Anthropic Cosmological Principle_ (Oxford University Press, Oxford, New York, 1986).
  * Aguirre (2001) A. Aguirre, Phys.Rev. D64, 083508 (2001), arXiv:astro-ph/0106143 [astro-ph] .
  * Harnik  _et al._ (2006) R. Harnik, G. D. Kribs, and G. Perez, Phys.Rev. D74, 035006 (2006), arXiv:hep-ph/0604027 [hep-ph] .
  * Adams (2008) F. C. Adams, JCAP 0808, 010 (2008), arXiv:0807.3697 [astro-ph] .
  * Note (2) An interesting additional note is that, here Hogan cites Hoyle, F., D. N. F. Dunbar, W. A. Wenzel, and W. Whaling, 1953, Phys. Rev. 92, 1095. This has been cited several times in different papers, sometimes with Phys. Rev. Lett. as the source, and a few times with the title “A state in C12 predicted from astrophysical evidence”. However, we failed to find any such article and would appreciate any information on this reference. The nearest match was D. N. F. Dunbar, R. E. Pixley, W. A. Wenzel, and W. Whaling, 1953, Phys. Rev. 92, 649, an article on the resonance in 12C often dubbed the ‘Hoyle resonance’.
  * Hanhart  _et al._ (2008) C. Hanhart, J. R. Pelaez, and G. Rios, Phys.Rev.Lett. 100, 152001 (2008), arXiv:0801.2871 [hep-ph] .
  * Pelaez and Rios (2010) J. R. Pelaez and G. Rios, Phys.Rev. D82, 114002 (2010), arXiv:1010.6008 [hep-ph] .
  * Bedaque  _et al._ (2011) P. F. Bedaque, T. Luu, and L. Platter, Phys.Rev. C83, 045803 (2011), arXiv:1012.3840 [nucl-th] .
  * Berengut  _et al._ (2013) J. Berengut, E. Epelbaum, V. Flambaum, C. Hanhart, U.-G. Meissner, _et al._ , Phys.Rev. D87, 085018 (2013), arXiv:1301.1738 [nucl-th] .