# Time Matters: Scaling Laws for Any Budget

Itay Inbar
Google DeepMind
<EMAIL_ADDRESS>
&Luke Sernau
Google DeepMind
<EMAIL_ADDRESS>

###### Abstract

A primary cost driver for training large models is wall-clock training time.
We show that popular time estimates based on FLOPs are poor estimates, and
construct a more accurate proxy based on memory copies. We show that with some
simple accounting, we can estimate the training speed of a transformer model
from its hyperparameters. Combined with a scaling law curve like Chinchilla,
this lets us estimate the final loss of the model. We fit our estimate to real
data with a linear regression, and apply the result to rewrite Chinchilla in
terms of a model’s estimated training time as opposed to the amount of
training data. This gives an expression for the loss in terms of the model’s
hyperparameters alone. We show that this expression is accurate across a wide
range of model hyperparameter values, enabling us to analytically make
architectural decisions and train models more efficiently.

##  1 Introduction

The final quality of a language model is constrained by the number of
parameters and the amount of data it was trained on. Remarkably, these two
parameters alone are often sufficient to estimate the final performance of the
model. Kaplan et al. [1] explored this phenomenon, predicting that loss curves
during pretraining could be written as a linear combination of a term
dependent on the number of the parameters and one dependent on the dataset
size. Hoffmann et al. [2] refined this estimate, improving the estimation of
the coefficients and introducing a bias term to capture the inherent
perplexity of language.

While these estimates are useful for large-scale models, small and mid-sized
models are not at risk of running out of pretraining data. Instead, the
limiting factor is the cost of training, a figure which is primarily driven by
a model’s size and speed. This suggests that instead of trading off model size
and dataset size, we should be trading off architectural hyperparameters
within the model that affect its throughput. On a fixed budget, a faster model
will be able to see more tokens than a slow one.

In this work we assume a fixed training time, and ask what hyperparameters we
should pick to maximize the final performance of the model. We start by
estimating the throughput of the model (tokens per second) in terms of the
number of FLOPs and memory copies, both of which can be directly calculated
from the model’s hyperparameters.

Kaplan et al. [1] mentions a parameterization in terms of compute
requirements, but their estimate is based on FLOPs, which we will show are a
weak predictor of runtime. Instead, we show that memory copies are a much
stronger predictor. This predictor, while simplistic, is powerful enough to
accurately predict the loss in terms of the hyperparameters of the model.

This new framing lets us estimate the final loss of a model without training
it, given only the model hyperparameters and the desired training time. We
show that this method produces accurate predictions across a wide range of
hyperparameter values, and makes useful predictions about which
hyperparameters should be used in order to maximize training efficiency.

We evaluate our findings over 1,535 different decoder-only transformer models
configurations ranging from $300K$ to $310M$ parameters and trained over the
C4 dataset [3]. We achieve an $r^{2}$ of 0.9 when predicting their final loss
using our refined scaling law. This is the same $r^{2}$ we get when using the
traditional Chinchilla scaling law. In other words, we are able to estimate
the final loss with the same accuracy whether we use Chinchilla scaling laws
on empirical runtimes or simply estimate them from hyperparameters.

##  2 The parameter equivalence principle

The core intuition motivating this work is the observation that large models
are not particularly sensitive to their hyperparameters, provided we hold the
total parameter count constant. This idea was discussed in [1], but due to its
importance we capture it in a form of an equivalence principle.

######  Principle 1 (The Parameter Equivalence Principle).

Above a certain scale, the final loss of a transformer is primarily a function
of how many parameters there are, not where they are in the model.

One straightforward implication is that we ought to be able to predict the
final loss using only parameter count and number of training tokens, as
earlier scaling laws did. But another often overlooked implication is that
models of the same size that allocate their parameters differently compete
primarily on speed.

If it is not feasible for one model to have vastly lower loss than another via
architectural improvements, we should instead choose architectures that
optimize for training speed, allowing them to consume as many tokens during
training time as it can. We can show this in practice by means of a scaling
law.

##  3 Estimating linear scaling law Coefficients

The original scaling law from Hoffmann et al. [2], predicts the final training
loss of a language model in terms of its parameters count $N$ and the number
of tokens it was trained upon $D$.

$L(N,D)=\frac{A}{N}^{\alpha}+\frac{B}{D}^{\beta}+E$ (1)

Hoffmann et al. [2] and Besiroglu et al. [4] derive different coefficient
values with the most extreme difference being their linear data coefficient
$B$. Rather than enter into this debate, we simply take the exponents from
[2], and fit our own linear coefficients $A$, $B$ and $E$ using linear
regression on the model loss. This was done by iterating over 1,535 different
decoder-only transformer models’ hyperparmeters configurations trained from
scratch for three hours each on the C4 Dodge et al. [3] dataset. Models sizes
vary from $319K$ to $310M$ parameters. Note that we constrained our
experiments to models that can be trained on a single TPU to avoid confounders
from inter-chip communication. We experimented with model hyperparameters of
embed sizes ranging from $2^{5}$ to $2^{10}$, number of layers ranging from
$3$ to $8$, MLP width ranging from $2^{8}$ to $2^{14}$, number of attention
heads ranging from $2^{1}$ to $2^{7}$, and a fixed vocabulary size of 8,000.
We trained on a mesh of 4 hosts x 8 chips/host of TPU V5 chips with no model
sharding.

Figure 1: Loss predictions over different trained models via Chinchilla, using
our linear coefficients.

The results show a very good fit ($r^{2}=0.9$), with coefficients $A=195.76$,
$B=182.52$, and $E=2.34$. Note that these are different from the values quoted
in either paper [2][4].

Using the values from the papers presents a very different story. The below
table compares the scaling law fitting measurements on our data using the
different papers coefficients, with our computed coefficients serving as a
baseline.

Table 1: Comparing scaling coefficients Fitted slope | Fitted intercept | $\mathbf{r^{2}}$ | Source
---|---|---|---
0.48 | 1.50 | 0.9 | Hoffmann et al.
0.43 | 1.70 | 0.9 | Besiroglu et al.
0.48 | 1.57 | 0.9 | Besiroglu et al. table params
1.00 | 0.00 | 0.9 | Hoffmann et al. with our A,B,E

Both Chinchilla papers underestimate the loss by a factor of more than two on
our data ($\text{slope}<0.5$). We take this as evidence that these
coefficients are perhaps highly sensitive to the details of the setup,
possibly explaining the discrepancy between the papers. Nonetheless, we were
able to achieve very good fit with a linear rescaling of their predictions
($r^{2}=0.9$ in all cases), suggesting the exponents are more robust.

##  4 Equations for estimating the speed of a model

In order to use scaling laws to estimate the loss of a model, we need to know
how big the model is and how much data will it be able to train over. The
former is a straightforward exercise in accounting, but the latter is more
nuanced.

We fix the amount of time we have to train the model to some constant $T$,
which is measured in seconds. If we can estimate how long each training step
takes for a given model, we can work out how many tokens (data) will be
processed by that model in time $T$.

It is tempting to imagine that we could estimate the model training speed just
by adding up the number of FLOPS. But as we will show, the runtime of the
model is actually driven by data copying, not the actual computation.

The amount of data copying depends on a wide variety of factors, from the
hardware to the architecture to the compiler. We do not attempt to account for
all of these factors here but take as a simplifying assumption that every
matrix multiplication requires a copy proportional to the size of its
operands.

Specifically, for a standard (as defined in [5]) decoder-only transformer
architecture, we derived equations for the number of parameters a model has
(PARAMS), the number of memory loads the model will need to make in a single
pass (MEMCPYS), and the number of operations the model will do in a single
pass (FLOPS):

$\displaystyle\text{PARAMS}\left(d,n,v,w\right)=vd+nd\left(8+2w+4d\right)+nw$
$\displaystyle\text{MEMCPYS}\left(d,n,s,v,w\right)=2vd+2sv+ns\left(w+2hs\right)+2nd\left(w+4s+2d\right)$
$\displaystyle\text{FLOPS}\left(d,n,s,v,w\right)=2svd+2dns\left(w+2d+s\right)+nhs^{2}$

Where the parameters are defined as:

$\qquad d=\text{embedding dimension}$

$\qquad n=\text{number of layers}$

$\qquad s=\text{sequence length}$

$\qquad v=\text{vocabulary size}$

$\qquad w=\text{MLP width}$

$\qquad h=\text{number of heads}$

The full details of the derivation of the above equations can be found in
appendix A.

Using a linear combination of the above equations we can now compute the total
number of seconds per training step (TIME) as:

$\text{TIME}\left(d,n,s,v,w\right)=c_{1}\text{MEMCPYS}\left(d,n,s,v,w\right)+c_{2}\text{FLOPS}\left(d,n,s,v,w\right)+c_{3}$
(2)

Where $c_{1}$, $c_{2}$, and $c_{3}$ are coefficients determined by linear
regression (see Section 5). Note that dividing the total number of seconds per
training step (i.e., TIME) by the number of seconds we are training upon
(i.e., $T$) would yield the total number of training steps (i.e., D in (1)).

Finally, we can estimate the total loss by plugging the above term into the
Chinchilla scaling law in order to derive an equation dependent on model
training speed.

$\hat{L}\left(d,n,s,v,w\right)=E+\frac{A}{\text{PARAMS}\left(d,n,v,w\right)^{\alpha}}+B\left(\frac{\text{TIME}\left(d,n,s,v,w\right)}{T}\right)^{\beta}$
(3)

Following our findings in Section 3, we take the original $\alpha$ and $\beta$
as in [2] and use our own fitted linear coefficients for $A$, $B$ and $E$.

##  5 Estimating the throughput

The throughput of a model is defined to be $1/\text{TIME}$, where TIME is
defined in Equation (2). In order to fully specify equation (2) we need to
determine its linear coefficients $c_{1}$, $c_{2}$, and $c_{3}$.

We conduct a large scale ($N=3,556$) sweep over model hyperparameters trained
for 5 minutes, just long enough to accurately determine the number of tokens
per second they process. We applied linear regression over the data to
determine $c_{1}$, $c_{2}$, and $c_{3}$.

We trained models of sizes varying from $277K$ parameters to $972M$
parameters. We experimented with model hyperparameters of embed sizes ranging
from $2^{5}$ to $2^{12}$, number of layers ranging from $1$ to $8$, MLP width
ranging from $2^{8}$ to $2^{15}$, number of heads ranging from $2^{0}$ to
$2^{7}$, and a fixed vocabulary size of 8,000. As in previous experiments we
trained on a mesh of 4x8 TPU V5 chips with no model sharding.

Figure 2: Estimating runtime with equation (2)

The results show an overall $r^{2}$ of $0.74$, with a much tighter fit for
slower (i.e. bigger) models. For fast models, confounding factors like
compiler optimizations start to matter, affecting the quality of the fit.

It is worth evaluating the importance of the different terms(i.e. FLOPS,
MEMCPY) in Equation (2). Previous work by both [2] and [4] utilized only the
FLOPS counting to derive their scaling laws. We show that MEMCPY is a stronger
predictor, and can account for essentially all of the explanatory power on its
own.

Figure 3: Runtime prediction ablation

##  6 Putting it all together

We now have an equation that estimates the number of tokens that the model
will consume from its hyperparameters. We also have an exact expression for
the number of parameters in such a model, PARAMS. Our tuned Chinchilla
equation relates these two quantities to estimate the final loss (3). In
figure 4, we show the results of this estimation, applied to the data from
Section 3.

Figure 4: Predicted vs actual loss. The prediction is made using only the
hyperparameters.

Notice that the graph is largely indistinguishable from Figure 1, including
the quality of fit $r^{2}=0.9$. While there is some error in the Chinchilla
equation’s predictions, there is essentially no additional error from using
our estimates in place of the empirical values.

(a) Empirical data consumption

(b) Estimated

Figure 5: Chinchilla using empirical data consumption vs estimated (ours)

##  7 Better loss with faster models

Figure 6: Negative gradient of the loss, projected such that the arrows point
in directions of constant parameter count. That is, following an arrow
corresponds to moving in the direction of steepest decent subject to the
constraint that the parameter count is held constant.

We can use these equations to make specific predictions about how we should
size our models. Figure 6 shows the negative gradient of the loss with respect
to each of our hyperparameters, projected to be along level curves of the
parameter count. Following each arrow brings you to another model with the
same parameter count but a lower predicted loss.

We can see that increasing the embed size at the expense of the other
hyperparameters is favorable throughout the plotted region. hidden_size is
particularly disincentivized. This suggests that we should take our MLPs to be
narrow, and our models to be somewhat shallow, in exchange for much larger
embed size.

##  8 Conclusion

Understanding what hyperparameters lead to the strongest model performance is
a vital part of model design. We’ve shown that the final loss of a model can
be accurately predicted by turning the question on its head. Instead of asking
for the most data efficient hyperparameters, we simply ask which
hyperparameters make the model the fastest. This leads to a new scaling law
based on hyperparameters alone. In the long run, the faster model will tend to
win.

We demonstrated this effect across a wide variety of model sizes, and showed
that we can accurately predict the model’s loss from its hyperparameters,
simply by estimating how many memory copies will take place. Crucially, this
is a stronger predictor than approaches based on FLOPs. What’s more, it allows
us to make specific predictions about which hyperparameters to use during
model design.

However, we do not consider the effects of model sharding, or the effects of
scale beyond a few hundred million parameters. We regard these as fruitful
areas of exploration for future work.

## References

  * Kaplan et al. [2020] Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei.  Scaling laws for neural language models.  _arXiv preprint arXiv:2001.08361_ , 2020.
  * Hoffmann et al. [2022] Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al.  Training compute-optimal large language models.  _arXiv preprint arXiv:2203.15556_ , 2022.
  * Dodge et al. [2021] Jesse Dodge, Maarten Sap, Ana Marasović, William Agnew, Gabriel Ilharco, Dirk Groeneveld, Margaret Mitchell, and Matt Gardner.  Documenting large webtext corpora: A case study on the colossal clean crawled corpus.  _arXiv preprint arXiv:2104.08758_ , 2021.
  * Besiroglu et al. [2024] Tamay Besiroglu, Ege Erdil, Matthew Barnett, and Josh You.  Chinchilla scaling: A replication attempt.  _arXiv preprint arXiv:2404.10102_ , 2024.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.

##  Appendix A Equations

###  A.1 FLOPS derivation

In order to compute the total number of FLOPS in our transformers decoders
stack we begin by counting the FLOPS needed for each step in a transformer
block.

Table 2: Transformer FLOPS $3sd^{2}$ | Produce $Q$, $K$ and $V$ inside the attention
---|---
$h\left(s^{2}d/h\right)$ | Compute $QK^{T}$
$h\left(s^{2}\right)$ | Apply softmax
$h\left(s^{2}d/h\right)$ | Multiply by $V$
$sd^{2}$ | Projection layer to recombine the heads
$sdw+swd$ | The two layers of the MLP

We sum all of these components and multiply by the number of transformer
layers in our transformer stack.

The final piece of the puzzle is to add the embedding of the input and the
output. Both of which require $svd$ FLOPS. We add all of these terms together
and simplify.

###  A.2 MEMCPYS derivation

We begin by adding up the total amount of data being copied in a single
transformer block. We approximate the number of memory copies needed for each
matmul as the size of the input matrices for each operation.

Table 3: MEMCPYS $3\left(sd+d^{2}\right)$ | Produce $Q$, $K$ and $V$ inside the attention
---|---
$h\left(sd/h+sd/h\right)$ | Compute $QK^{T}$
$h\left(s^{2}\right)$ | Apply softmax
$h\left(s^{2}+sd/h\right)$ | Multiply by $V$
$sd+d^{2}$ | Projection layer to recombine the heads
$sd+dw+sw+wd$ | MLP

Again we sum all of the above components and multiply by the number of
transformer layers in our transformer stack.

Finally, we add the embedding of the input and the output, both of which
require $v*d+s*v$ memory copies. Summing all of this together and simplifying
yields our MEMCPY equation.

###  A.3 PARAMS derivation

We begin with the per-layer parameters. We note the extra vector term to
account for the bias term accompanying each matrix.

Table 4: PARAMS $3\left(d+d^{2}\right)$ |  $Q$, $K$ and $V$ inside the attention
---|---
$4d$ | Layer norms before and after the attention
$d^{2}+d$ | Projection layer to recombine the heads
$2dw+d+w$ | MLP matrices

In similar fashion, we left with summing all of the above components and
multiplying by the number of transformer layers in our transformer stack.
The final piece of the puzzle is to add a $2d$ vector for the norm layer after
the transformers as well as the embedding matrix used for embedding the input
and the output, matrix of size $vd$. Summing all of these and simplifying
yields our PARAMS equation.


# One nine availability of a Photonic Quantum Computer on the Cloud toward HPC
integration

Nicolas Maring1, Andreas Fyrillas12, Mathias Pont1, Edouard Ivanov1,
Eric Bertasi1, Mario Valdivia1 and Jean Senellart1  1Quandela, Massy, France
Email<EMAIL_ADDRESS>2Centre for Nanosciences and
Nanotechnologies, CNRS, Université Paris-Saclay, Palaiseau, France

###### Abstract

The integration of Quantum Computers (QC) within High-Performance Computing
(HPC) environments holds significant promise for solving real-world problems
by leveraging the strengths of both computational paradigms. However, the
integration of a complex QC platform in an HPC infrastructure poses several
challenges, such as operation stability in non-laboratory like environments,
and scarce access for maintenance. Currently, only a few fully-assembled QCs
currently exist worldwide, employing highly heterogeneous and cutting-edge
technologies. These platforms are mostly used for research purposes, and often
bear closer resemblance to laboratory assemblies rather than production-ready,
stable, and consistently-performing turnkey machines. Moreover, public cloud
services with access to such quantum computers are scarce and their
availability is generally limited to few days per week. In November 2022, we
introduced the first cloud-accessible general-purpose quantum computer based
on single photons. One of the key objectives was to maintain the platform’s
availability as high as possible while anticipating seamless compatibility
with HPC hosting environment. In this article, we describe the design and
implementation of our cloud-accessible quantum computing platform, and
demonstrate one nine availability (92%) for external users during a six-month
period, higher than most online services. This work lay the foundation for
advancing quantum computing accessibility and usability in hybrid HPC-QC
infrastructures.

###### Index Terms:

Quantum Computing, Hosting, Availability

##  I Introduction

The first quantum computer available on the cloud was introduced by IBM in May
2016 [1]. IBM’s quantum computing platform, known as IBM Quantum Experience,
provided researchers and developers with remote access to a small-scale
quantum computer through the cloud. This marked an important milestone in the
democratization of quantum computing, allowing users to experiment and explore
quantum algorithms without the need for physical access to the hardware [2].
Since then, several other companies and organizations have also made quantum
computers accessible through direct cloud-based platforms or third-party
services [3], contributing to the growth and accessibility of quantum
computing.

However, these resources have relatively low availability, with system
uptime111as evident from the device availability dashboards of different
platforms, for instance AWS Console AWS Braket > Device (https://us-
east-1.console.aws.amazon.com/braket/home#/devices), or Microsoft Azure
AzureQuantum > Providers generally ranging from one day to four days per week,
which falls far below the high availability standards set for classical
resources on the same platforms (e.g., three to five nines). For all the
providers, the reasons for this limited availability are likely rooted in
multiple challenges:

  * •

Limited hardware resources: quantum computers are in the early stages of their
development, and the number of available platforms is limited. This scarcity
leads to a high demand for access to quantum computers, resulting in limited
availability for users.

  * •

Complex hardware maintenance: quantum computers require specialized
maintenance and calibration procedures to ensure their proper functioning.
These processes are time-consuming and require dedicated, highly trained
teams.

  * •

Error rates and stability: quantum computers are highly sensitive to
environmental noise and susceptible to errors. Achieving high stability and
low error rates over a long period of time is a significant challenge in
quantum computing. When errors occur, it requires troubleshooting and
recalibration, resulting in periods of unavailability.

  * •

Research and development focus: many quantum computers available through cloud
services are still in the research and development phase (TRL 4 or above).
These systems may undergo frequent updates, modifications, and
experimentation, leading to intermittent periods of unavailability as new
features or improvements are implemented.

We released in November 2022 the first single-photon based computing platform
named Ascella [4], and provided free online access through the Quandela Cloud
service [5]. In this work, we present the architecture of Ascella and the
specific challenges to making the platform openly accessible. Inspired from
the “early days” of classical computers [6], we describe the protocols,
processes, and adaptations implemented to achieve a first nine availability of
the service. Furthermore, we provide an overview of the remaining challenges
expected for scaling up to larger systems and a roadmap towards achieving
high-performance computing (HPC) standards. Lastly, we emphasize the necessity
of transparency in availability and real-time performance metrics for the
possible and efficient integration of QC into HPC environment.

##  II Architecture and key operation challenges

Ascella’s hardware, as shown in Fig. 1.a, is based on an on-demand high-
brightness source of pure and indistinguishable single photons [7]. This
source, made from a semiconductor quantum dot (behaving as “an artificial
atom”) embedded in a micropillar cavity, requires both an applied bias voltage
and a finely tuned excitation laser to emit pure single photons [8]. These are
then interfaced with a deterministic optical demultiplexer [9], generating and
initializing up to 6 path-encoded (dual-rail) qubits. Photon arrival times are
synchronized using fiber delay lines, ensuring on-chip photons wave packets
overlap. The qubits are manipulated in a 12-mode fully re-configurable
universal photonic integrated circuit (UPIC) [10], and measured with high
efficiency single-photon detectors and a post-processing unit.

To operate the machine, tasks are sent by the user to Quandela Cloud through
standard REST API using the python-based open-source framework Perceval [11].
The user can either specify an available high-level primitive (for instance
photonic Variational Quantum Eigensolver on a specific molecule), a photonic
circuit, a gate-based circuit or a unitary transformation to be applied to a
specified input state containing $1$ to $6$ photonic qubits, and optional
post-selection criteria. Output photon coincidences are then acquired up to
the desired sample number and data sample results are sent back to the user to
be interpreted in the context of the running algorithm.

Figure 1:  Architecture of the Ascella QPU a. A quantum-dot single-photon
source (SPS) device at 5 is operated at 80 repetition rate. An active optical
demultiplexer (DMX) followed by fibered delays converts the train of single
photons into 6 photons arriving simultaneously to the 12-mode universal
photonic integrated chip (UPIC). Photons are detected at the chip output by
superconducting nanowire single-photon detectors (SNSPD) and detection times
are processed by a correlator. A full software stack controls the unitary
matrix implemented on the chip through the voltages applied on 126 thermo-
optic phase shifters and the photonic input state according to the job
requested. b. The setup may be integrated into a server-compatible structure.
c. Continuous benchmark over 120 days of operation. CNOT is the average of a
simple CNOT fidelity run at different position on the photonic chip, and RB36
operates a full tomography of a CNOT gate through 36 sampling taks as detailed
in [4] ” _Benchmarking logic gates_ ”. d. Perceval users may send jobs to the
Quandela Cloud consisting in either a high-level primitive, like a variational
quantum eigensolver for a molecule, or an input state along with a photonic
circuits (PC), a gate-based circuit (GB), or a unitary matrix (U). The job is
first processed by a CPU, which computes the optimized voltages for the on-
chip phase shifters from our compilation and transpilation process. Finally,
the user receives the collected single-photon and coincidence counts after the
computation on the quantum processing unit (QPU).

As common in a research facility, the full set-up was designed and assembled
on an optical table, as a mosaic of in-house and off-the-shelf components
(laser, cryostats, mirrors, lenses, optical fibers, etc…). Active components
are connected to a single desktop computer through USB and controlled through
heterogeneous proprietary software stacks and ad-hoc communication protocols.
Most of the components are operated at room temperature with standard air
conditioning. To enable a stable operation of the cloud-available system, it
is necessary to consider hardware faults and fluctuating performance. Indeed,
the quantum computing system consists of various components susceptible to
faults such as communication errors, feed-back control instability, or
internal software bugs. The main challenge is that these failures may not
break the entire system flow, leading to the production of inconsistent
results which the end-user may not even be aware of. Anticipating hardware
faults is essential to ensuring the stability and reliability of the system.
Furthermore, hardware components involved in generating and manipulating
qubits are sensitive to environmental fluctuations. For instance, spikes in
current can modify the voltage bias applied to the single-photon source,
typical room temperature fluctuations can alter the refractive index of delay
fibers, causing slight delays between photons and polarization rotations, and
room humidity fluctuations can impact free-space modules. Managing these
fluctuations and their impact on the system’s performance is crucial for
maintaining consistent and reliable results.

The seamless user experience of our cloud-accessible computer leads to
additional constraints. For a cloud-based service, users outside our close
collaborators must be taken into account: it is imperative to maintain
consistent availability and responsiveness, even during non-working hours, as
well as programmed and openly accessible system updates and maintenance
windows. Ensuring seamless communication, setting realistic expectations, and
providing prompt support is vital to user satisfaction.

##  III Processes

To enhance hardware reliability, expedite maintenance response time, and
elevate user experience, we implemented a series of processes.

###  III-A Monitoring and alerts

Every hour, Ascella undergoes a few minutes automated calibration where 92
system parameters and metrics, such as optical transmissions, single-photon
purity and indistinguishability [12], and modules temperatures, are measured.
Note that this calibration period is counted as down-time in our availability
estimation. The metrics are sent to an aggregation server called Minotaur,
which includes a MQTT server Mosquitto and the standardized data aggregation
server Prometheus. The logs are exported to a log data
platform222https://www.ovhcloud.com/en/logs-data-platform/ utilizing the
graylog standard. All metrics are visualized on a dashboard using a Grafana
visualization server and have been displayed in break rooms and high-
visibility locations. Based on predetermined thresholds on key metrics, the
Grafana alert engine sends anomaly notifications to the cloud operations team.
This tool plays a crucial role in rapidly diagnosing failures and facilitating
swift maintenance actions.

Finally, we collect “usage data”: all the technical information related to
user tasks sent to the platform (in average $32\,000$ per month in the first 6
months). Usage data is stored preserving only timestamp, time of operation,
compiled circuit, input photon configuration and task numeric results. All
these logs are kept with a 6 month retention period and used for system
stability and predictive maintenance.

###  III-B Methodological reaction and feedback control system

Our approach to mitigating the fluctuating performance discussed in Sec. II –
and minimizing the frequency of raised alerts – is based on an agile-like
incident management methodology. This methodology relies on incident triage
and tracking, down-time transparency, swift reaction, root cause analysis,
monitoring (see Sec. III-A), automated testing and continuous improvement. In
particular, we use active feedback loops on a number of parameters in the
opto-electronic system. These loops are kept as simple as possible to ensure
robustness against measurement noise: they mostly consist of coordinate
descents with fixed-step line searches. To proactively anticipate performance
decline and expedite corrective measures, we have identified a set of five
pivotal parameters in the qubit generation and manipulation modules (see Fig.
1.a) which we actively provide feedback on; to maximise brightness (i.e. the
photon extraction rate), we control: (1) the single-photon source bias
voltage, (2) the photon polarization at the input of the optical
demultiplexer, and (3) the photon polarization at each of the six inputs of
the Ascella QPU; brightness is furthermore stabilized by (4) maintaining the
excitation laser power at a fixed, pre-determined value; finally, to maintain
photon arrival-time synchronicity (which is strongly related to two-qubit gate
fidelity), (5) the fiber delays at the output of the optical demultiplexer are
kept at constant temperature in a thermalized box. Each anomaly exceeding the
scope of this stabilisation is used to develop more advanced simulation
modules for each hardware component of the system, thereby increasing the
automated test coverage.

Any set-up improvements and upgrades to the QPU hardware or software stack are
then reviewed and discussed on a weekly basis, and implemented during the
weekly 4-hour maintenance period.

While photonic integrated chips guarantee unprecedented compactness and
stability for optical setups, they are subject to imperfections: passive
phases [13], phase shifter thermal crosstalk [14], asymmetric directional
couplers [15] and inhomogeneous input/output losses [16]. We developed a
machine learning-based characterization method that is able to measure these
imperfections, which we can then effectively compensate. Our custom
compilation process optimizes the applied on-chip phase shifts to enhance the
fidelity of the implemented operation acting on the qubits. Subsequently, our
tailored transpilation process translates the phase shifts into voltages used
to control the photonic chip (see Fig. 1.a), taking crosstalk in particular
into account. The characterization process can take up to 24 hours for the
Ascella UPIC and must be executed every time a new single-photon source is
installed, but it is also compatible with on-the fly updates, which is being
implemented, leveraging “usage data” acquired during its operation time.

###  III-C Continuous benchmark and output analysis

In addition to monitoring physical system parameters as presented in Section
III-A, we also introduced very early high-level benchmarks in order to catch
unanticipated system failures, and monitor the evolution of global performance
over long periods of time. These benchmarks are designed to non-trivially use
the full capacity of the processor to avoid blind spots in the monitoring: we
compute 1-, 2- and 3-qubit gate fidelity using the protocol described in [4]
through respectively 4, 36 and 340 independent measurements. To avoid delay of
user tasks, these measurements are launched by subset of measurements ran
every 2 hours using period with no user activity. The resulting fidelity
therefore aggregates measurements over long periods, incidentally highlighting
any time-dependent fluctuations. Figure 1.c shows for instance the monitoring
of CNOT [17] through 12 full measurements per day over 120 days.

Finally, we exploit “usage data” by comparing processor output with
simulators333The possibility to fully simulate output of QPUs is currently an
opportunity for system understanding and development that we have to leverage
fully since it will disappear when number of qubits will grow: for such
photonic quantum computer, we expect to reach this point where output will not
be anymore simulable for about 20 photons. also available on the platform and
using this comparison to progressively update noise models of the simulator
and again to detect new types of system failures.

###  III-D User communication

Up-to-date information regarding the hardware’s performance is directly
available on Quandela Cloud: the current status (available, in deployment, in
calibration, in maintenance, or unreachable) of each computing platform, the
number of tasks currently queuing, a detailed description of the hardware are
available online, and the latest characterization of the key performance
index. The performance metrics include the indistinguishability of the
generated photons (Hong–Ou–Mandel visibility), the total transmittance of the
setup and the single-photon purity ($\mathcal{P}=1-g^{(2)}(0)$). The physical
constraints imposed by the hardware are also provided to the user. All this
information can be used by the end-user to analyze outcome of his own
algorithm. For example, the Hong–Ou–Mandel measurement is expected to have a
direct impact on CNOT fidelity [18]. For each request, the user has therefore
the ability to use this imperfection information for getting a better
interpretation of his results.

Scheduled and unplanned maintenance periods are communicated through a News
tab.

Finally, a support link is provided in the Cloud interface and dedicated
personnel respond to any incoming request regardless of the user subscription
plan.

##  IV Conclusion

This 6-month-long operational experience has been highly beneficial and
instructive: it has allowed us to identify numerous points of failure and
implement automated procedures, significantly reducing the need for “human in
the loop” for troubleshooting and manual calibration. We have collected
stability and performance benchmarks throughout the entire period, which has
been invaluable for training better models and predicting complex failures or
system erosion. Throughout this process, we have maintained close contact with
end-users through direct and indirect communication.

Reaching the milestone of achieving first nine availability is important to
us, although it may be considered modest in comparison to HPC standards.
Nonetheless, we believe that it should become a minimal requirement. Reaching
the second nine remains out of reach for the moment, given the inherent
complexity of the hardware and the challenges in applying redundancy, which is
a key aspect of achieving the highest availability in classical architectures.

In addition to operational improvements, we have also made significant
progress in developing a data center-compatible solution adhering to 19-inch
rack standards, as depicted in Figure 1.b. This integrated system incorporates
an alignment-free and mechanically vibration-insensitive pigtailed single-
photon source, along with single photon detectors, within a compact and
industrial 4 cryostat. The system architecture primarily relies on in-house
opto-electronic modules, which utilizes a unified and centralized
communication protocol and power management system. This design facilitates
easy remote maintenance, incorporates auto-alignment and optimization
functions, and eliminates the need for specialized operator intervention.
Monitoring is an integral part of the system’s design rather than an
additional layer. The system has been intentionally built to incorporate
monitoring capabilities from the ground up, ensuring efficient performance and
reliability. Its compact and modular integration significantly mitigates
environmental instabilities, allows for easy parts replacement, and serves as
the foundation for future scaled-up versions of quantum computers.

While the journey towards full integration of Quantum Computers in full HPC
environment and achieving quantum advantage in real-world applications just
started and is expected to open new challenges, we advocate for a pragmatic
incremental bottom-up approach. By transparently sharing each step with end-
users, we can collectively acquire the necessary knowledge and methodologies
for operating reliable quantum processors. High availability of the quantum
processing unit (QPU) is crucial for successful integration into HPC
environments but is not sufficient. Simultaneously, this integration process
should be strongly guided by application-specific benchmarks to evaluate
performance, optimize algorithms, and drive further advancements in quantum
computing tailored to specific use cases.

## References

  * [1] IBM, “The quantum experience. url http://www. research.ibm.com/quantum/.”
  * [2] S. J. Devitt, “Performing quantum computing experiments in the cloud,” _Physical Review A_ , vol. 94, no. 3, p. 032329, 2016.
  * [3] V. Hassija, V. Chamola, V. Saxena, V. Chanana, P. Parashari, S. Mumtaz, and M. Guizani, “Present landscape of quantum computing,” _IET Quantum Communication_ , vol. 1, no. 2, pp. 42–48, 2020.
  * [4] N. Maring, A. Fyrillas, M. Pont, E. Ivanov, P. Stepanov, N. Margaria, W. Hease, A. Pishchagin, T. H. Au, S. Boissier _et al._ , “A general-purpose single-photon-based quantum computing platform,” _arXiv preprint arXiv:2306.00874_ , 2023.
  * [5] Quandela, “Quandela cloud,” 2022, https://cloud.quandela.com.
  * [6] J. Gray and D. P. Siewiorek, “High-availability computer systems,” _Computer_ , vol. 24, no. 9, pp. 39–48, 1991.
  * [7] N. Somaschi, V. Giesz, L. De Santis, J. Loredo, M. P. Almeida, G. Hornecker, S. L. Portalupi, T. Grange, C. Anton, J. Demory _et al._ , “Near-optimal single-photon sources in the solid state,” _Nature Photonics_ , vol. 10, no. 5, pp. 340–345, 2016.
  * [8] S. E. Thomas, M. Billard, N. Coste, S. C. Wein, Priya, H. Ollivier, O. Krebs, L. Tazaïrt, A. Harouri, A. Lemaitre, I. Sagnes, C. Anton, L. Lanco, N. Somaschi, J. C. Loredo, and P. Senellart, “Bright Polarized Single-Photon Source Based on a Linear Dipole,” _Physical Review Letters_ , vol. 126, no. 23, jun 2021.
  * [9] M. Pont, R. Albiero, S. E. Thomas, N. Spagnolo, F. Ceccarelli, G. Corrielli, A. Brieussel, N. Somaschi, H. Huet, A. Harouri _et al._ , “Quantifying n-photon indistinguishability with a cyclic integrated interferometer,” _Physical Review X_ , vol. 12, no. 3, p. 031033, 2022.
  * [10] C. Taballione, R. van der Meer, H. J. Snijders, P. Hooijschuur, J. P. Epping, M. de Goede, B. Kassenberg, P. Venderbosch, C. Toebes, H. v. d. Vlekkert _et al._ , “A 12-mode universal photonic processor for quantum information processing,” _arXiv preprint arXiv:2012.05673_ , 2020.
  * [11] N. Heurtel, A. Fyrillas, G. de Gliniasty, R. Le Bihan, S. Malherbe, M. Pailhas, E. Bertasi, B. Bourdoncle, P.-E. Emeriau, R. Mezher _et al._ , “Perceval: A software platform for discrete variable photonic quantum computing,” _Quantum_ , vol. 7, p. 931, 2023.
  * [12] H. Ollivier, S. Thomas, S. Wein, I. M. de Buy Wenniger, N. Coste, J. Loredo, N. Somaschi, A. Harouri, A. Lemaitre, I. Sagnes _et al._ , “Hong-ou-mandel interference with imperfect single photon sources,” _Physical Review Letters_ , vol. 126, no. 6, p. 063602, 2021.
  * [13] Y. Yang, Y. Ma, H. Guan, Y. Liu, S. Danziger, S. Ocheltree, K. Bergman, T. Baehr-Jones, and M. Hochberg, “Phase coherence length in silicon photonic platform,” _Optics Express_ , vol. 23, no. 13, p. 16890–16902, Jun 2015.
  * [14] M. Milanizadeh, D. Aguiar, A. Melloni, and F. Morichetti, “Canceling thermal cross-talk effects in photonic integrated circuits,” _Journal of Lightwave Technology_ , vol. 37, no. 4, p. 1325–1332, Feb 2019.
  * [15] S. Bandyopadhyay, R. Hamerly, R. Hamerly, and D. Englund, “Hardware error correction for programmable photonics,” _Optica_ , vol. 8, p. 1247–1255, Oct 2021.
  * [16] M. Li, C. Zou, G. Guo, and X. Ren, “Effect of unbalanced and common losses in quantum photonic integrated circuits,” vol. 15, p. 092701, Sep 2017.
  * [17] T. C. Ralph, N. K. Langford, T. Bell, and A. White, “Linear optical controlled-not gate in the coincidence basis,” _Physical Review A_ , vol. 65, no. 6, p. 062324, 2002.
  * [18] O. Gazzano, M. Almeida, A. Nowak, S. Portalupi, A. Lemaître, I. Sagnes, A. White, and P. Senellart, “Entangling quantum-logic gate operated with an ultrabright semiconductor single-photon source,” _Physical review letters_ , vol. 110, no. 25, p. 250501, 2013.


* [36] Z. Wei, W. Yuan, S. Li, J. Yuan, and D. W. K. Ng, “Transmitter and receiver window designs for orthogonal time-frequency space modulation,” IEEE Trans. Commun., vol. 69, no. 4, pp. 2207-2223, Apr. 2021.
  * [37] S. Coleri, M. Ergen, A. Puri and A. Bahai, “Channel estimation techniques based on pilot arrangement in OFDM systems,” _IEEE Tran. Broadcas._ , vol. 48, no. 3, pp. 223-229, Sept. 2002.
  * [38] E. Björnson and L. Sanguinetti, “Making cell-free massive MIMO competitive with MMSE processing and centralized implementation,” _IEEE Trans. Wireless Commun.,_ vol. 19, no. 1, pp. 77-90, Jan. 2020.
  * [39] Z. Wang, J. Zhang, B. Ai, C. Yuen and M. Debbah, “Uplink performance of cell-free massive MIMO with multi-antenna users over jointly-correlated Rayleigh fading channels,” _IEEE Trans. Wireless Commun.,_ Accepted, 2022.
  * [40] Guidelines for Evaluation of Radio Interface Technologies for IMTAdvanced, document Rep. ITU 2135-1, 2009.
  * [41] F. Zhang, Matrix Theory: Basic Results and Techniques. Springer Science & Business Media, 2011.


the expression for $\delta_{j}=e^{\beta\Delta_{j}}$ is also associated with
the toy-model, but since $\delta_{j}$ is defined as a ratio of binding
affinities, a purely exponential representation of $\delta_{j}$ is arguably
general. Taking $\gamma_{j}=(\beta E_{V})^{3/2}e^{-\beta E_{j}}$ was primarily
for analytic convenience and real systems would have a different thermal
dependence for $\gamma_{j}$.. Using the grid-assembly expression for
$\gamma_{j}$ amounts to using the second equation in Eq.(89) to define the
critical temperature. The Jensen-inequality based derivation that shows how
this equation limits $N$ is similar to that which leads to Eq.(94).
Ultimately, we find

$N<(\beta_{\text{crit}}E_{V})^{3/2}e^{\beta_{\text{crit}}(\bar{\Delta}+\bar{E})},$
(95)

where $\bar{E}\equiv\sum_{j=1}^{D}n_{j}E_{j}/N$. Thus the implication is the
same as that which follows Eq.(94): Ligand-receptor systems which have a large
number of ligands and where optimal binding is physically important, should
also have large average binding energies for those ligands.

Eq.(86) gives the thermal condition under which fully optimal ligand-receptor
binding can be achieved, but there is not a universal reason for why such a
bound configuration would be functionally optimal in all biomolecular systems.
Perhaps the microstate where all ligands are in their optimal receptor sites
is too rigid to be biophysically useful, and thus it is preferred if the
system exists in a partially-bound state where most (but not all) of the bound
ligands are in their optimal configuration. Thus the value in this formalism
might not exist in its ability to predict the temperature at which the system
is in a fully bound state, but rather in its power to predict and categorize
binding properties above this temperature. For this value, the system-type
distinction is the major biophysical utility of this formalism: Ligand-
receptor systems are inherently combinatorial, but it seems that if nature
were to evolve such systems so as to avoid the kinetic traps of suboptimal
bindings, it would engineer ligands to have sufficiently high binding
affinities to their optimal-receptors as to easily achieve the combinatorics-
limit Eq.(87) at in vivo temperatures.

## 8 Discussion

We began this work with the goal of using combinatorics to model how the
competition between distinct ligands affects the ligands’ equilibrium binding
properties to receptors. Before developing a physical model of such binding,
we needed to solve a modified version of a well known problem on derangements.
After solving this problem and using it to compute a general partition
function, we were able to derive implicit formulas for the average number of
bound ligands and the average number of optimally bound ligands for each type
(Eq.(79)). Narrowing our focus to the case of matched receptor and ligand
populations (i.e., $r_{j}=n_{j}$ for all $j$), we were able to derive the
condition under which all ligands are bound to their optimal receptors
(Eq.(86)).

The observables defined in Eq.(79) allowed us to compute the equilibrium
binding properties of our ligand-receptor system at any temperature (provided
we have thermal dependences for $\gamma_{j}$ and $\delta_{j}$ for all $j$).
But the model also provided softer biophysical predictions. Eq.(94) and
Eq.(95) can be seen as such predictions, labeled as ”soft” because rather than
explicitly predicting that some observable has a value, they predict that, in
order for a system condition to be satisfied, an observable cannot exceed a
certain value. The main implication of these inequalities is that the larger
the binding energy advantage for optimal contacts, the more particles the
system can have and still be capable of achieving the fully optimally bound
state. Thus, systems with more ligands should also have larger optimal binding
energies supposing such bindings have functional importance. Note that Eq.(95)
was derived for the simple model where free-ligands were taken to be point
particles, but an analogous inequality could be obtained for more physically
informed values of $\gamma_{j}$.

For the task of quantitatively modeling optimal ligand-receptor binding, there
are many limitations to the introduced model.

First, we are only partly modeling the combinatorics of our system, a
combinatorics that can be characterized by derangements from a pre-defined
”correct” or optimal configuration. A more general and flexible model would
not a priori define such a configuration and would instead have an interaction
matrix defining how various ligand types interact with various receptor types,
and would then use this matrix to determine the optimal matchings between
ligands and receptors. However, it is not clear how such a more general
formalism would be soluble. We started Sec. 3 by representing our problem in
terms of a general interaction matrix (in Eq.(13)), but then had to make the
simplifying assumption of a pre-defined optimal configuration in order to put
the partition function in an analytically tractable form. If we wanted to
maintain analytical tractability while pursuing a more general model, we would
likely need a similar simplifying assumption.

Second, to derive our various optimal binding conditions (Eq.(57), Eq.(72),
and Eq.(82)), we assumed that the system contained the same number of receptor
sites as ligands for each type (i.e., $n_{j}=r_{j}$ for all $j$). This
assumption was convenient for the framework of derangements, but does not
match a realistic biophysical scenario for which the number of ligands and
receptors are not necessarily equal. Thus, a reasonable extension to the
Biophysical Implications Sec. 7 would be to find the analytical conditions
that define optimal binding for the case of unmatched receptor and ligand
populations of each type. However, the general form of the observables in this
system (given in Eq.(79)) do not make the assumption of $n_{j}=r_{j}$ for all
$j$, and thus can still be used to numerically compute the number of bound and
unbound-particles of various types as a function of temperature.

Third, we did not discuss how the spatial organization of receptors can affect
the propensity of ligands to bind to them. Instead we assumed that all
receptors were on equal footing in terms of spatial accessibility, and all
binding variances could be encoded into the collection of parameters
$\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{D})$ and
$\boldsymbol{\delta}=(\delta_{1},\ldots,\delta_{D})$. How spatial organization
affects binding could likely be better modeled through a matrix-based
interaction scheme where spatially occluded receptors have reduced binding
affinities, but, as previously mentioned, such an interaction scheme is less
analytically tractable than that assumed in the paper.

The first and the third of these realistic additions would complicate the
combinatorial simplicity of what we have presented here, so much so that an
entirely new model would likely have to be erected in its place. As in all
modeling, building such a new model would require us to balance the increase
in relevance from an incorporation of more realistic properties with a
concomitant loss in transparency from the new model’s relative insolubility.

Finally, the model we developed does not answer the implicit question posed in
the introduction. In Sec. 1, we first motivated our initial steps towards the
general model presented in Sec. 3 by noting that many ligand-receptor systems
exhibit many-to-many interactions in which each type of ligand can bind to
many receptors and vice versa. This fact led us to build a model in which even
when ligands interacted optimally with a few receptors, these ligands still
had the potential to interact with sub-optimal receptors. However, the results
of Sec. 6.1 suggest that if optimal ligand-receptor binding (where each ligand
binds to its energetically optimal receptor) is functionally important in the
system, then that binding has to be highly specific to combat its
combinatorial disadvantage. So although our model was motivated by multi-
specific ligand binding, the model’s results primarily pertain to highly-
specific ligand binding and the thermal conditions needed to achieve such
binding. Therefore left over from the introduction is the question of how does
one quantitatively model the ligand-receptor binding underlying the
combinatorial codes that make multi-specific systems so necessary in
pharmacology and real biomolecular environments. Building such a model would
likely require us to make use of an interaction matrix since it is only
through such a matrix that one can encode different sets of ligands
interacting optimally with different sets of receptors.

As an ending remark, we note that although Sec. 4.2 explored a limiting case
to our more general model, the model introduced there can also stand alone as
a probabilistic model of derangements. For a state space of permutations of a
list with repeated elements (e.g., all of the ways to permute the characters
and spaces in a sentence), we assume that there is a single permutation where
the elements are said to be in their ”correct order.” Let the parameter
$w_{i}$ be proportional to the probability that an element of type $i$ is in
its correct position. Then, the partition function we previously computed in
Eq.(59) becomes the weighted sum of states

$X_{\boldsymbol{n}}(\boldsymbol{w})=\int^{\infty}_{0}dx\,e^{-x}\,\prod_{k=1}^{D}\left(w_{k}-1\right)^{n_{k}}\,L^{(0)}_{n_{k}}\left(\frac{x}{1-w_{k}}\right),$
(96)

where we took $n_{j}=r_{j}$ since we are considering permutations of a list
where there are matched populations of correct positions and elements. Eq.(96)
is essential to the computation of two quantities: First, the probability that
$\ell_{i}$ elements of type $i$, for $i=1,\ldots,D$, are in their correct
positions; Second, the average number of elements that are in their correct
positions. Respectively, these quantities are

$P_{\boldsymbol{\ell}}=\frac{G_{\boldsymbol{n-\ell}}}{X_{\boldsymbol{n}}(\boldsymbol{w})}\prod_{k=1}^{D}\binom{n_{k}}{\ell_{k}}w_{k}^{\ell_{k}},\qquad\langle\ell_{\text{tot}}\rangle=\sum_{j=1}^{D}n_{j}w_{j}\frac{X_{\boldsymbol{n}_{j}}(\boldsymbol{w})}{X_{\boldsymbol{n}}(\boldsymbol{w})},$
(97)

where $G_{\boldsymbol{n}}$ (defined in Eq.(2) and related to
$B_{\boldsymbol{r},\boldsymbol{k}}$ through
$G_{\boldsymbol{n}}=B_{\boldsymbol{n},\boldsymbol{n}}$) is the generalized
derangement formula, and $\boldsymbol{n}_{j}$ is $\boldsymbol{n}$ with 1
subtracted from the $j$th component:
$\boldsymbol{n}_{j}=(n_{1},\ldots,n_{j}-1,\ldots,n_{D})$. Aside from deriving
them directly, we can obtain the expressions in Eq.(97) by translating the
Boltzmann factor expressions in Sec. 4.2 into the language of probability
weights by taking $\delta_{j}\to w_{j}$. Analogous to Eq.(93), using Jensen’s
inequality, we find that a necessary, but not sufficient, condition for the
average $\langle\ell_{\text{tot}}\rangle$ to be equal to $N$ (i.e., for all
elements to be in their correct order) is

$1\leq\sum_{j=1}^{D}\frac{n_{j}\ln w_{j}}{N\ln N}.$ (98)

This probabilistic model of derangements would be relevant in a computing
context where one is interested in the probability of various derangements
that differ from a given sequence by a fixed number of elements.

## 9 Conclusion

We introduced a derangement model of ligand-receptor interactions. The model
advances the subject of ligand-receptor modeling in two directions: First, it
concretely frames the question of how a finite collection of ligands can
compete for a finite collection of receptors. Such a question is important
because ligands in real biomolecular systems always exist in crowded
environments with other types of ligands, and, with the multi-specificity of
ligand-receptor binding, such crowded environments ordinarily exhibit ligands
struggling to bind to their optimal sites in a sea of suboptimal ones. Second,
to better organize the counting of the system’s microstates, the model
introduces a combinatorial problem (and its solution) as the foundational
framework for such systems. This approach to building biomolecular models by
beginning with combinatorial questions appears to be generalizable to similar
contexts since biomolecular systems often contain finite numbers of particles
where particles can only interact in precise ways that can often be defined by
some combinatorial set.

## 10 Acknowledgements

The author thanks Michael Brenner and Rostam Razban for helpful comments on
the work.

## 11 Author’s Statements

  * •

Financial Support: This research received no specific grant from any funding
agency, commercial or nonprofit sectors.

  * •

Conflict of Interests Statement: The author has no conflicts of interest to
disclose.

  * •

Ethics Statement: This research did not require ethical approval.

## 12 Supplementary Code

The code used to generate all figures is found in the repository
https://github.com/mowillia/LigandReceptor.

##  Appendix A Simulations

To simulate our system and create the plots in Fig. 6(b), Fig. 7(b), Fig.
8(b), and Fig. 9 we implemented a Metropolis Hastings algorithm [Kra06]. The
code for recreating these figures is linked to in Sec. 12. Here we review the
salient parts of the implementation.

For these simulations, we needed to define a microstate, the probability of
transitions between microstates, and the types of transitions between
microstates.

###  A.1 Microstate Definition

A microstate of our system was defined by two lists: one representing the
collection of unbound particles, and the other representing particles bound to
their various binding sites. The particles themselves were denoted by unique
strings and came in multiple copies according to the system parameters. For
example, a system with $D=3$ types of particles with $n_{1}=2$, $n_{2}=3$, and
$n_{3}=1$ could have a microstate defined by
$\texttt{unbound\\_particles}=[A_{2},A_{2},A_{3}]$ and
$\texttt{bound\\_particles}=[A_{1},-,A_{2},-,A_{1},-]$ where ”$-$” in the
bound list stands for an empty binding site.

Since the number of optimally bound particles was an important observable for
the system, we also needed to define the optimal binding configuration for the
microstates. Such an optimal configuration was chosen at the start of the
simulation and was defined as a microstate with no unbound particles and all
the bound particles in a particular order. For example, using the previous
example, we might define the optimal binding configuration as
$\texttt{optimal\\_bound\\_config}=[A_{1},A_{1},A_{2},A_{2},A_{2},A_{3}]$, in
which case the number of optimally bound particles of each type in
$\texttt{bound\\_particles}=[A_{1},-,A_{2},-,A_{1},-]$ is $m_{1}=1$,
$m_{2}=1$, and $m_{3}=0$. The number of bound particles of each type is
$k_{1}=2$, $k_{2}=1$, and $k_{3}=0$. We note that the order of the elements in
unbound_particles is not physically important, but, since the number of
optimally bound particles is an important observable, the order of the
elements in bound_particles is physically important.

For these simulations, the temperature-normalized energy of a microstate with
$k_{i}$ bound particles of type $i$ and $m_{i}$ optimally bound particles of
type $i$ was defined as

$\beta
E(\boldsymbol{k},\boldsymbol{m})=\sum_{i=1}^{D}(m_{i}\ln\delta_{i}+k_{i}\ln\gamma_{i}),$
(99)

where $\boldsymbol{k}=(k_{1},k_{2},\ldots,,k_{D})$,
$\boldsymbol{m}=(m_{1},m_{2},\ldots,m_{D})$, $\gamma_{i}$ is the binding
affinity, and and $\delta_{i}$ is the optimal-binding affinity.

###  A.2 Transition Probability

For transitioning between microstates, we allowed for three different
transition types: Particle binding to a site; particle unbinding from a site;
permutation of two particles in two different binding sites. Particle binding
and unbinding both occur in real physical systems, but permutation of particle
positions is unphysical. This latter transition type was included to ensure an
efficient-in-time sampling of the state space. For simulations of equilibrium
systems it is valid to include physically unrealistic transition types as long
as the associated transition probabilities obey detailed balance.

At each time step, we first randomly selected one of the three transition
types with equal probability for each type, then randomly selected the final
proposed microstate given the initial microstate, and finally computed the
probability that said proposal was accepted. By the Metropolis Hastings
algorithm [Kra06], the probability that the transition is accepted is given by

$\text{acceptance
prob}(\text{init}\to\text{fin})=\min\left\\{1,e^{-\beta(E_{\text{fin}}-E_{\text{init}})}\frac{\pi(\text{fin}\to\text{init})}{\pi(\text{init}\to\text{fin})}\right\\},$
(100)

where $\beta$ is inverse temperature, $E_{\text{init}}$ is the energy of the
initial microstate state, and $E_{\text{fin}}$ is the energy of the final
microstate, with energy defined in Eq.(99). The quantity
$\pi(\text{init}\to\text{fin})$ is the probability of randomly proposing the
final microstate state given the initial microstate state and
$\pi(\text{fin}\to\text{init})$ is defined similarly. The ratio
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})$ varied for each
transition type.

###  A.3 Transition Types and Examples

Below we give examples of the three types of transitions along with the value
of the ratio $\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})$ in
each case. In the following, $N_{f}$ and $N_{b}$ represent the number of free
particles and the number of bound particles, respectively, before the
transition.

  * •

Particle Binding to Site: One particle was randomly chosen from the
unbound_particles list and placed in a randomly chosen empty site in the
bound_particles list.
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})=(N_{b}+1)^{-1}/(N_{f}^{-1}\times
N_{f}^{-1})={N_{f}^{2}}/{(N_{b}+1)}$.
Example:
$\texttt{unbound\\_particles}=[A_{2},A_{2},A_{3}]$ and
$\texttt{bound\\_particles}=[A_{1},-,A_{2},-,A_{1},-]$
(transitions to) $\\\ \texttt{unbound\\_particles}=[A_{2},A_{3}]$ and
$\texttt{bound\\_particles}=[A_{1},A_{2},A_{2},-,A_{1},-]$;
Transition weight:
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})=9/4$

  * •

Particle Unbinding from Site: One particle was randomly chosen from the
bound_particles list and placed in the unbound_particles list.
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})=((N_{f}+1)^{-1}\times(N_{f}+1)^{-1})/(N_{b}^{-1})={N_{b}}/{(N_{f}+1)^{2}}$.
Example:
$\texttt{unbound\\_particles}=[A_{2},A_{2},A_{3}]$ and
$\texttt{bound\\_particles}=[A_{1},-,A_{2},-,A_{1},-]$
(transitions to) $\\\ \texttt{unbound\\_particles}=[A_{2},A_{2},A_{3},A_{2}]$
and $\texttt{bound\\_particles}=[A_{1},-,-,-,A_{1},-]$;
Transition weight:
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})=3/16$

  * •

Particle Permutation: Two particles were randomly selected from the
bound_particles list, and their positions in the list were switched.
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})=1.$
Example:
$\texttt{unbound\\_particles}=[A_{2},A_{2},A_{3}]$ and
$\texttt{bound\\_particles}=[A_{1},-,A_{2},-,A_{1},-]$
(transitions to) $\\\ \texttt{unbound\\_particles}=[A_{2},A_{2},A_{3}]$ and
$\texttt{bound\\_particles}=[A_{2},-,A_{1},-,A_{1},-]$;
Transition weight:
$\pi(\text{fin}\to\text{init})/\pi(\text{init}\to\text{fin})=1$

For impossible transitions (e.g., particle binding when there are no free
particles) the probability for accepting the transition was set to zero. At
each temperature, the simulation was run for anywhere from 10,000 to 30,000
time steps depending on observed convergence, and the final $2\%$ of the time
steps were used to compute ensemble averages of $\langle k\rangle$ and
$\langle m\rangle$. These simulations were repeated five times, and each point
in Fig. 6(b), Fig. 7(b), Fig. 8(b), and Fig. 9 represents the average $\langle
k\rangle$ and $\langle m\rangle$ over these five runs.

##  Appendix B Consistency Checks for $B_{\boldsymbol{r},\boldsymbol{k}}$

In this section, we affirm the various consistency checks for
$B_{\boldsymbol{r},\boldsymbol{k}}$. The first check ensures that
$B_{\boldsymbol{r},\boldsymbol{k}}$ has the proper limiting case when
$\boldsymbol{r}$ has components with value $1$. The second check ensures that
$B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}$ has the
correct normalization when summed over all possible values of
$\boldsymbol{m}$.

###  B.1 Checking Eq.(8)

We want to check that

$\sum_{k_{1}=0}^{1}\ldots\sum_{k_{D}=0}^{1}B_{\boldsymbol{r}_{0},\boldsymbol{k}}\delta(K,k_{1}+\cdots+k_{D})=\sum_{j=0}^{D}(-1)^{j}\binom{D}{j}\binom{D-j}{K-j}^{2}(K-j)!,$
(101)

where $\boldsymbol{r}_{0}=(1,\ldots,1)$. We start from Eq.(4) (reproduced here
for convenience):

$\displaystyle B_{\boldsymbol{r},\boldsymbol{k}}$
$\displaystyle=\sum_{j_{1}=0}^{k_{1}}\cdots\sum_{j_{D}=0}^{k_{D}}\binom{r_{1}}{j_{1}}\cdots\binom{r_{D}}{j_{D}}(-1)^{j_{1}+\cdots+j_{D}}\binom{r_{1}-j_{1}+\cdots+r_{D}-j_{D}}{k_{1}-j_{1}+\cdots+k_{D}-j_{D}}$
$\displaystyle\qquad\qquad\times\frac{(k_{1}-j_{1}+\cdots+k_{D}-j_{D})!}{(k_{1}-j_{1})!\cdots(k_{D}-j_{D})!}.$
(102)

Given $j_{i}\leq k_{i}\leq r_{i}$, if $r_{i}=1$ then we must have
$(k_{i}-j_{i})!=1$. We also have $\binom{r_{i}}{j_{i}}=\binom{1}{j_{i}}=1$.
Thus Eq.(4) becomes

$\displaystyle B_{\boldsymbol{r}_{0},\boldsymbol{k}}$
$\displaystyle=\sum_{j_{1}=0}^{k_{1}}\cdots\sum_{j_{D}=0}^{k_{D}}(-1)^{j_{1}+\cdots+j_{D}}\binom{D-j_{1}-\cdots-
j_{D}}{K-j_{1}-\cdots-j_{D}}$
$\displaystyle\qquad\qquad\times{(K-j_{1}-\cdots-j_{D})!},$ (103)

where we defined $K\equiv\sum_{i=1}^{D}k_{i}$. With the identity,

$1=\sum_{J=0}^{K}\delta(J,j_{1}+\cdots+j_{D}),$ (104)

we can introduce a Kronecker delta and obtain

$\displaystyle B_{\boldsymbol{r}_{0},\boldsymbol{k}}$
$\displaystyle=\sum_{J=0}^{K}\sum_{j_{1}=0}^{k_{1}}\cdots\sum_{j_{D}=0}^{k_{D}}(-1)^{j_{1}+\cdots+j_{D}}\binom{D-j_{1}-\cdots-
j_{D}}{K-j_{1}-\cdots-j_{D}}$
$\displaystyle\qquad\qquad\times{(K-j_{1}-\cdots-
j_{D})!}\,\delta(J,j_{1}+\cdots+j_{D})$
$\displaystyle=\sum_{J=0}^{K}(-1)^{J}\binom{D-J}{K-J}(K-J)!\sum_{j_{1}=0}^{k_{1}}\cdots\sum_{j_{D}=0}^{k_{D}}\delta(J,j_{1}+\cdots+j_{D}).$
(105)

Isolating the summation over the Kronecker delta yields

$\displaystyle\sum_{j_{1}=0}^{k_{1}}\cdots\sum_{j_{D}=0}^{k_{D}}\delta(J,j_{1}+\cdots+j_{D})$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dz}{z^{J+1}}\prod_{i=1}^{D}\sum_{j_{i}=0}^{k_{i}}z^{j_{i}}$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dz}{z^{J+1}}(1+z)^{k_{1}+\cdots+k_{D}}$
$\displaystyle=\binom{K}{J}.$ (106)

Thus, we obtain

$\displaystyle B_{\boldsymbol{r}_{0},\boldsymbol{k}}$
$\displaystyle=\sum_{J=0}^{K}(-1)^{J}\binom{D-J}{K-J}(K-J)!\binom{K}{J}.$
(107)

Performing the final summation in Eq.(101), we obtain

$\displaystyle\sum_{k_{1}=0}^{1}\ldots\sum_{k_{D}=0}^{1}B_{\boldsymbol{r}_{0},\boldsymbol{k}}\delta(K,k_{1}+\cdots+k_{D})$
$\displaystyle=\sum_{J=0}^{K}(-1)^{J}\binom{D-J}{K-J}(K-J)!\binom{K}{J}$
$\displaystyle\qquad\qquad\times\sum_{k_{1}=0}^{1}\ldots\sum_{k_{D}=0}^{1}\delta(K,k_{1}+\cdots+k_{D})$
$\displaystyle=\sum_{J=0}^{K}(-1)^{J}\binom{D-J}{K-J}(K-J)!\binom{K}{J}\binom{D}{K},$
(108)

and with the identity

$\binom{K}{J}\binom{D}{K}=\binom{D}{J}\binom{D-J}{K-J},$ (109)

we have

$\sum_{k_{1}=0}^{1}\ldots\sum_{k_{D}=0}^{1}B_{\boldsymbol{r}_{0},\boldsymbol{k}}\delta(K,k_{1}+\cdots+k_{D})=\sum_{J=0}^{K}(-1)^{J}\binom{D-J}{K-J}^{2}(K-J)!\binom{D}{J},\quad\checkmark$
(110)

as expected.

Why does this result make sense? When
$\boldsymbol{r}=(1,\ldots,1)\equiv\boldsymbol{r}_{0}$, the vector
$\boldsymbol{k}$ can only have elements of $1$ or $0$. Thus
$B_{\boldsymbol{r}_{0},\boldsymbol{k}}$ represents the number of ways to
completely derange a particular collection of objects, defined by
$\boldsymbol{k}$, out of $D$ unique objects. In order to find $b_{D,K}$, the
number of ways to completely derange $K$ objects, we need to count and sum the
number of derangements for all collections of objects. Therefore, to obtain
$b_{D,K}$ we need to sum $B_{\boldsymbol{r}_{0},\boldsymbol{k}}$ over all
possible values of $\boldsymbol{k}$ such that $\sum_{j}k_{j}=K$.

###  B.2 Checking equivalence between Eq.(10) and Eq.(11)

We want to find a reduced form for

$I_{\boldsymbol{r},\boldsymbol{k}}=\sum_{m_{1}=0}^{k_{1}}\cdots\sum_{m_{D}=0}^{k_{D}}\binom{r_{1}}{m_{1}}\cdots\binom{r_{D}}{m_{D}}B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}.$
(111)

The expression for $B_{\boldsymbol{r},\boldsymbol{k}}$ is

$B_{\boldsymbol{r},\boldsymbol{k}}=\frac{1}{(\sum_{i}\alpha_{i})!}\int^{\infty}_{0}dx\,e^{-x}\prod_{i=1}^{D}(-1)^{k_{i}}x^{\alpha_{i}}L_{k_{i}}^{(\alpha_{i})}(x),$
(112)

Noting that $\alpha_{i}=r_{i}-k_{i}=(r_{i}-m_{i})-(k_{i}-m_{i})$, from
Eq.(111) and Eq.(112), we find

$I_{\boldsymbol{r},\boldsymbol{k}}=\frac{1}{(\sum_{i}\alpha_{i})!}\int^{\infty}_{0}dx\,e^{-x}\,x^{\sum_{i}\alpha_{i}}\prod_{i=1}^{D}\sum_{m_{i}=0}^{k_{i}}\binom{r_{i}}{m_{i}}(-1)^{k_{i}-m_{i}}L_{k_{i}-m_{i}}^{(\alpha_{i})}(x).$
(113)

Next, we make the change of variables

$q_{i}\equiv k_{i}-m_{i}.$ (114)

We then have

$I_{\boldsymbol{r},\boldsymbol{k}}=\frac{1}{(\sum_{i}\alpha_{i})!}\int^{\infty}_{0}dx\,e^{-x}x^{\sum_{i}\alpha_{i}}\prod_{i=1}^{D}\sum_{q_{i}=0}^{k_{i}}\binom{k_{i}+\alpha_{i}}{k_{i}-q_{i}}(-1)^{q_{i}}L_{q_{i}}^{(\alpha_{i})}(x).$
(115)

To simplify this result we need to find an identity for

$\sum_{q=0}^{k}\binom{k+\alpha}{k-q}U^{q}L_{q}^{(\alpha)}(x).$ (116)

Expanding the generalized Laguerre polynomial $L^{(\alpha)}_{q}(x)$ according
to its definition, we obtain

$\displaystyle\sum_{q=0}^{k}\binom{k+\alpha}{k-q}U^{q}L_{q}^{(\alpha)}(x)$
$\displaystyle=\sum_{q=0}^{k}\binom{k+\alpha}{k-q}U^{q}\sum_{i=0}^{q}\binom{q+\alpha}{q-i}\frac{(-1)^{i}}{i!}x^{i}.$
(117)

More identity wrangling gives us

$\binom{k+\alpha}{k-q}\binom{q+\alpha}{q-i}=\binom{k+\alpha}{k-i}\binom{k-i}{q-i},$
(118)

and thus Eq.(117) becomes

$\displaystyle\sum_{q=0}^{k}\binom{k+\alpha}{k-q}U^{q}L_{q}^{(\alpha)}(x)$
$\displaystyle=\sum_{i=0}^{k}\binom{k+\alpha}{k-i}\frac{(-1)^{i}}{i!}x^{i}\sum_{q=i}^{k}U^{q}\binom{k-i}{q-i}$
$\displaystyle=\sum_{i=0}^{k}\binom{k+\alpha}{k-i}\frac{(-1)^{i}}{i!}x^{i}U^{i}\sum_{q=i}^{k}U^{q-i}\binom{k-i}{q-i}$
$\displaystyle=\sum_{i=0}^{k}\binom{k+\alpha}{k-i}\frac{(-1)^{i}}{i!}x^{i}U^{i}(1+U)^{k-i},$
(119)

which, by the definition of the Laguerre polynomial yields

$\sum_{q=0}^{k}\binom{k+\alpha}{k-q}U^{q}L_{q}^{(\alpha)}(x)=(1+U)^{k}L_{k}^{(\alpha)}\left(\frac{xU}{1+U}\right).$
(120)

Also, from the definition of the Laguerre polynomial, we can show

$\lim_{\lambda\to
0}\lambda^{k}L_{k}^{(\alpha)}\left(\frac{x}{\lambda}\right)=(-1)^{k}\frac{x^{k}}{k!}.$
(121)

Therefore, with Eq.(120) and Eq.(121), we have

$\sum_{q=0}^{k}\binom{k+\alpha}{k-q}(-1)^{q}L_{q}^{(\alpha)}(x)=\lim_{U\to-1}(1+U)^{k}L_{k}^{(\alpha)}\left(\frac{xU}{1+U}\right)=\frac{x^{k}}{k!}.$
(122)

Inserting this result into Eq.(115) gives us

$\displaystyle I_{\boldsymbol{r},\boldsymbol{k}}$
$\displaystyle=\frac{1}{(\sum_{i}\alpha_{i})!}\int^{\infty}_{0}dx\,e^{-x}x^{\sum_{i}\alpha_{i}}\frac{x^{k_{1}+\cdots+k_{D}}}{k_{1}!\cdots
k_{D}!}$ $\displaystyle=\frac{1}{(\alpha_{1}+\cdots+\alpha_{D})!k_{1}!\cdots
k_{D}!}\int^{\infty}_{0}dx\,e^{-x}x^{r_{1}+\cdots+r_{D}},$ (123)

which, with the Gamma function definition, yields

$I_{\boldsymbol{r},\boldsymbol{k}}=\frac{(r_{1}+\cdots+r_{D})!}{(r_{1}-k_{1}+\cdots+r_{D}-k_{D})!k_{1}!\cdots
k_{D}!}.$ (124)

Why does this result make sense? From one perspective, the quantity
$\binom{r_{1}}{m_{1}}\cdots\binom{r_{D}}{m_{D}}B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}$
is the number of ways to choose $m_{j}$ out of $r_{j}$ positions to contain
their correct elements while the remaining $k_{j}-m_{j}$ elements are
completely deranged with respect to their $r_{j}-m_{j}$ remaining correct
positions, for $j=1,\ldots,D$. If we sum over all possible values of $m_{j}$
we should obtain the number of ways to arrange (i.e., not only derange)
$k_{j}\leq r_{j}$ objects for $j=1,\ldots,D$ across a total of
$r_{1}+\ldots+r_{D}$ lattice sites.

On the other hand, if we are trying to calculate the number of ways to arrange
$k_{j}$ objects of type $j$ for $j=1,\ldots,D$ across $r_{1}+\cdots+r_{D}$
lattice sites, we can use the language of multinomials. Say that the objects
represent ”filled” sites on the lattice and the spaces between objects are
”empty” sites. Then there are $k_{1}+\cdots+k_{D}$ filled sites (of which
$k_{i}$ are identical for each $i$) and $r_{1}-k_{1}+\cdots+r_{D}-k_{D}$ empty
sites. Finding the number of ways to arrange the objects amongst the
$r_{1}+\cdots+r_{D}$ sites is equivalent to finding the total number of ways
to order this collection of filled and empty sites while correcting for
equivalent orderings due to reordering the positions of the same type of site.
Including both filled and empty sites there is a total of $r_{1}+\cdots+r_{D}$
sites amongst which we have $r_{1}-k_{1}+\cdots+r_{D}-k_{D}$ ”copies” of empty
sites, $k_{1}$ copies of filled sites of type 1, $k_{2}$ copies of filled
sites of type 2, $\ldots$ and $k_{D}$ copies of filled sites of type $D$.
Counting the number of ways to order the $r_{1}+\cdots+r_{D}$ sites and
correcting for the equivalent reorderings arising from the multiple copies of
various types of sites leads to the multionomial

$\frac{(r_{1}+\cdots+r_{D})!}{(r_{1}-k_{1}+\cdots+r_{D}-k_{D})!k_{1}!\cdots
k_{D}!},$ (125)

which we have shown is equivalent to the result written in terms of a
summation over $B_{\boldsymbol{r-m},\boldsymbol{k-m}}$.

##  Appendix C Deriving General Partition Function

In this section, we derive Eq.(24) from Eq.(22). First we note from the
definition of $B_{\boldsymbol{r},\boldsymbol{k}}$ in Eq.(6) that

$B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}=\frac{1}{(\sum_{j}\alpha_{j})!}\int^{\infty}_{0}dx\,e^{-x}x^{\sum_{j}\alpha_{j}}\prod_{j=1}^{D}(-1)^{k_{j}-m_{j}}L_{k_{j}-m_{j}}^{(\alpha_{j})}(x),$
(126)

where $\alpha_{j}=r_{j}-k_{j}$. Thus the summation over $\boldsymbol{m}$ in
Eq.(22) becomes

$\displaystyle\sum_{\boldsymbol{m}}B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}\prod_{j=1}^{D}\binom{r_{j}}{m_{j}}\delta_{j}^{m_{j}}$
$\displaystyle=\frac{1}{(\sum_{j}\alpha_{j})!}\int^{\infty}_{0}dx\,e^{-x}x^{\sum_{j}\alpha_{j}}$
$\displaystyle\qquad\times\prod_{j=1}^{D}\sum_{m_{j}=0}^{k_{j}}\binom{r_{j}}{m_{j}}\delta_{j}^{m_{j}}(-1)^{k_{j}-m_{j}}L_{k_{j}-m_{j}}^{(\alpha_{j})}(x),$
$\displaystyle=\frac{1}{(\sum_{j}\alpha_{j})!}\int^{\infty}_{0}dx\,e^{-x}x^{\sum_{j}\alpha_{j}}$
$\displaystyle\qquad\times\prod_{j=1}^{D}\sum_{q_{j}=0}^{k_{j}}\binom{k_{j}+\alpha_{j}}{k_{j}-q_{j}}\delta_{j}^{k_{j}-q_{j}}(-1)^{q_{j}}L_{q_{j}}^{(\alpha_{j})}(x),$
(127)

where we changed variables from $m_{j}$ to $q_{j}=k_{j}-m_{j}$ in the final
line. Using the following Laguerre polynomial identity (derived in Eq.(119))

$\sum_{q=0}^{k}\binom{k+\alpha}{k-q}U^{q}L_{q}^{(\alpha)}(X)=(1+U)^{k}L_{k}^{(\alpha)}\left(\frac{XU}{1+U}\right),$
(128)

we then have

$\sum_{m_{j}=0}^{k_{j}}\binom{k_{j}+\alpha_{j}}{k_{j}-q_{j}}\delta_{j}^{k_{j}-q_{j}}(-1)^{q_{j}}L_{q_{j}}^{(\alpha_{j})}(x)=(\delta_{j}-1)^{k_{j}}L_{k_{j}}^{(\alpha_{j})}\left(\frac{x}{1-\delta_{j}}\right).$
(129)

Therefore, the summation over $\boldsymbol{m}$ becomes

$\displaystyle\sum_{\boldsymbol{m}}B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}\prod_{j=1}^{D}\binom{r_{j}}{m_{j}}\delta_{j}^{m_{j}}$
$\displaystyle=\frac{1}{(\sum_{j}\alpha_{j})!}\int^{\infty}_{0}dx\,e^{-x}\prod_{j=1}^{D}x^{\alpha_{j}}(\delta_{j}-1)^{k_{j}}L_{k_{j}}^{(\alpha_{j})}\left(\frac{x}{1-\delta_{j}}\right)$
$\displaystyle\equiv C_{\boldsymbol{r},\boldsymbol{k}}$ (130)

For the summation over $\boldsymbol{k}$, we use the contour integral identity
for inverse factorial

$\frac{1}{(\sum_{j}\alpha_{j})!}=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z^{\sum_{j}\alpha_{j}+1}}e^{z},$ (131)

where $\Gamma$ is a closed contour about the origin, to obtain

$\displaystyle\sum_{\boldsymbol{k}}C_{\boldsymbol{r},\boldsymbol{k}}$
$\displaystyle\prod_{j=1}^{D}\frac{1}{(n_{j}-k_{j})!}\,\gamma_{j}^{k_{j}}=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z}e^{z}\int^{\infty}_{0}dx\,e^{-x}\prod_{j=1}^{D}\sum_{k_{j}=0}^{n_{j}}\left(\frac{x}{z}\right)^{\alpha_{j}}\frac{\gamma_{j}^{k_{j}}(\delta_{j}-1)^{k_{j}}}{(n_{j}-k_{j})!}L_{k_{j}}^{(\alpha_{j})}\left(\frac{x}{1-\delta_{j}}\right)$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z}e^{z}\int^{\infty}_{0}dx\,e^{-x}\prod_{j=1}^{D}\left(\frac{x}{z}\right)^{r_{j}}\sum_{k_{j}=0}^{n_{j}}\frac{1}{(n_{j}-k_{j})!}\left(\frac{z\gamma_{j}}{x}(\delta_{j}-1)\right)^{k_{j}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times
L_{k_{j}}^{(r_{j}-k_{j})}\left(\frac{x}{1-\delta_{j}}\right).$ (132)

Using the identity

$\sum_{k=0}^{n}\frac{Y^{k}}{(n-k)!}L_{k}^{(r-k)}(X)=Y^{n}\,L_{n}^{(r-n)}(X-Y^{-1}),$
(133)

proved in Appendix D, we find that the final partition function is

$\displaystyle{\cal Z}_{\boldsymbol{n},\boldsymbol{r}}$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z}\int^{\infty}_{0}dx\,e^{z-x}\prod_{j=1}^{D}(\gamma_{j}(\delta_{j}-1))^{n_{j}}\left(\frac{x}{z}\right)^{r_{j}-n_{j}}L_{n_{j}}^{(r_{j}-n_{j})}\left(\frac{x(z\gamma_{j}+1)}{z\gamma_{j}(1-\delta_{j})}\right)$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z}\int^{\infty}_{0}dx\,e^{z-x}\left(\frac{x}{z}\right)^{N_{R}-N_{L}}\prod_{j=1}^{D}(\gamma_{j}(\delta_{j}-1))^{n_{j}}L_{n_{j}}^{(r_{j}-n_{j})}\left(\frac{x(z\gamma_{j}+1)}{z\gamma_{j}(1-\delta_{j})}\right),$
(134)

where we defined $N_{R}\equiv\sum_{j=1}^{D}r_{j}$ and
$N_{L}\equiv\sum_{j=1}^{D}n_{j}$ as the total number of receptors and total
number of ligands, respectively, in the system. In deriving Eq.(134), we have
made no assumptions about the relative value of $r_{j}$ and $n_{j}$ for each
$j$ and therefore the result is valid for both $r_{j}<n_{j}$ and $r_{j}\geq
n_{j}$.

##  Appendix D Laguerre Polynomial Argument Summation Identity Eq.(133)

In this section, we prove the identity Eq.(133). Starting from the definition
of the generalized Laguerre polynomial Eq.(5), we have

$\displaystyle\sum_{k=0}^{n}\frac{Y^{k}}{(n-k)!}L_{k}^{(r-k)}(X)$
$\displaystyle=\sum_{k=0}^{n}\frac{Y^{k}}{(n-k)!}\sum_{j=0}^{k}\binom{k+r-k}{k-j}\frac{(-1)^{j}}{j!}X^{j}$
$\displaystyle=\sum_{k=0}^{n}\frac{Y^{k}}{(n-k)!}\sum_{\ell=0}^{k}\binom{r}{\ell}\frac{(-1)^{k-\ell}}{(k-\ell)!}X^{k-\ell}$
$\displaystyle=\sum_{\ell=0}^{n}\sum_{k=\ell}^{n}Y^{k}(-1)^{k-\ell}X^{k-\ell}\binom{r}{\ell}\frac{1}{(n-\ell)!}\binom{n-\ell}{k-\ell}$
(135)

where we changed summation variables ($\ell=k-j$) in the second line and
switched the order of summations in the third line. Separating the summations
and then using the binomial theorem, we obtain

$\displaystyle\sum_{k=0}^{n}\frac{Y^{k}}{(n-k)!}L_{k}^{(r-k)}(X)$
$\displaystyle=\sum_{\ell=0}^{n}\binom{r}{\ell}\frac{1}{(n-\ell)!}Y^{\ell}\sum_{k=\ell}^{n}(-1)^{k-\ell}(XY)^{k-\ell}\binom{n-\ell}{k-\ell}$
$\displaystyle=\sum_{\ell=0}^{n}\binom{r}{\ell}\frac{1}{(n-\ell)!}Y^{\ell}(1-XY)^{n-\ell}$
$\displaystyle=Y^{n}\sum_{\ell=0}^{n}\binom{r}{\ell}\frac{(-1)^{n-\ell}}{(n-\ell)!}(X-Y^{-1})^{n-\ell}$
$\displaystyle=Y^{n}\sum_{j=0}^{n}\binom{n+r-n}{n-j}\frac{(-1)^{j}}{j!}(X-Y^{-1})^{j}$
$\displaystyle=Y^{n}\,L^{(r-n)}_{n}(X-Y^{-1})$ (136)

where we changed variables to $j=n-\ell$ in the fourth line and used the
definition Eq.(5) in the final line.

##  Appendix E Gendered Dimer Assembly Equilibrium Conditions

Here we derive the equilibrium conditions Eq.(39) and Eq.(40) from the large
$N$ approximation conditions Eq.(36) and the observable definitions Eq.(37)
and Eq.(38). We start with the conditions defining $\bar{x}$ and $\bar{z}$:

$1=\sum_{j=1}^{D}\frac{\bar{z}\gamma_{j}+1}{\bar{x}+\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})},\qquad\bar{z}=\bar{x}\sum_{j=1}^{D}\frac{1}{\bar{x}+\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})}.$
(137)

The expressions for $\langle k_{j}\rangle$ and $\langle m_{j}\rangle$ are

$\displaystyle\langle k_{j}\rangle$
$\displaystyle=\frac{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})}{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})+\bar{x}}$
(138) $\displaystyle\langle m_{j}\rangle$
$\displaystyle=\frac{\bar{z}\gamma_{j}\delta_{j}}{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})+\bar{x}}.$
(139)

Next, we seek to eliminate the $\bar{z}$ and $\bar{x}$ dependence from these
observables. From Eq.(137), Eq.(138), and Eq.(139), we can show

$\bar{z}=N-\langle k\rangle,\qquad\bar{x}=N-\sum_{j=1}^{D}\langle
m_{j}\rangle(1-\delta_{j}^{-1}),$ (140)

where we used $\langle k\rangle=\sum_{j=1}^{D}\langle k_{j}\rangle$, and the
total particle number $N$ is equal to the number of particle types $D$ when
there is one copy per particle. With Eq.(138) and Eq.(139), we can also show

$\langle k_{j}\rangle-\langle
m_{j}\rangle(1-\delta_{j}^{-1}\rangle=\frac{\bar{z}\gamma_{j}\bar{x}}{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})+\bar{x}},$
(141)

and with the second equation in Eq.(137) and the first equation in Eq.(140),
we obtain

$\sum_{j=1}^{D}\frac{1}{\gamma_{j}}\Big{(}\langle k_{j}\rangle-\langle
m_{j}\rangle(1-\delta_{j}^{-1})\Big{)}=\Big{(}N-\langle k\rangle\Big{)}^{2},$
(142)

which is our first equilibrium condition.

From both equations in Eq.(137), we have

$\sum_{j=1}^{D}\frac{\bar{z}\gamma_{j}}{\bar{x}+\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})}=1-\sum_{j=1}^{D}\frac{1}{\bar{x}+\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})}=1-\bar{z}/\bar{x}.$
(143)

Using Eq.(143) in Eq.(141) and Eq.(139), we find, respectively

$\displaystyle\sum_{j=1}^{D}\Big{(}\langle k_{j}\rangle-\langle
m_{j}\rangle(1-\delta_{j}^{-1}\rangle\Big{)}$ $\displaystyle=\bar{x}-\bar{z}$
(144) $\displaystyle\sum_{j=1}^{D}\langle m_{j}\rangle\delta_{j}^{-1}$
$\displaystyle=1-\bar{z}/\bar{x}.$ (145)

With Eq.(144) and Eq.(145) and the second equation in Eq.(140), we then obtain

$\displaystyle\sum_{j=1}^{D}\langle m_{j}\rangle\delta_{j}^{-1}=\frac{\langle
k\rangle-\langle m\rangle+\sum_{j=1}^{D}\langle
m_{j}\rangle\delta_{j}^{-1}}{N-\langle m\rangle+\sum_{j=1}^{D}\langle
m_{j}\rangle\delta_{j}^{-1}},$ (146)

which is our second equilibrium condition.

##  Appendix F Derivation of Eq.(82)

In this section, we derive Eq.(82), the general thermal condition for the
fully optimally bound state. This state is defined as

$\langle k\rangle=\langle m\rangle=N-1+O(\delta^{-1}).$ (147)

where $N\equiv\sum_{j=1}^{D}n_{j}=\sum_{j=1}^{D}r_{j}$. To derive the desired
thermal condition, we will need some of our previous results. In particular,
we need the average number of optimally bound ligands of type $j$

$\langle m_{j}\rangle=\frac{n_{j}\delta_{j}}{\delta_{j}-1}\frac{\displaystyle
L_{n_{j}-1}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}\left(\bar{\phi}_{j}\right)},$ (148)

and the equations relating $\bar{z}$ and $\bar{x}$ to our observables:

$\bar{z}=N-\langle k\rangle,\qquad\bar{x}=N-\sum_{i=1}^{D}\langle
m_{i}\rangle(1-\delta_{i}^{-1}).$ (149)

From Eq.(147), we can infer that for the fully optimally bound state we have

$\langle m_{j}\rangle=n_{j}-O(D^{-1})+O(\delta_{j}^{-1}),$ (150)

where we are explicitly referencing the $O(D^{-1})$ terms in order to satisfy
$\sum_{j=1}^{D}O(D^{-1})=1$ required of Eq.(147). Using Eq.(150) in the second
equation in Eq.(149) yields

$\bar{x}=1+\sum_{j=1}^{D}n_{j}\delta^{-1}_{j}+O(\delta^{-1}),$ (151)

where we subsumed terms of order $O(D^{-1}\delta^{-1})$ into $O(\delta^{-1})$.

For the fully optimally bound configuration, ligands need to greatly favor
their optimal receptors in order to bind only to such receptors. Thus for the
fully optimally bound state, we can assume that the optimal-binding affinity
for each particle-type is much greater than 1: $\delta_{j}\gg 1$. From this
assumption we can take $\bar{\phi}_{j}$ defined as

$\bar{\phi}_{j}=\frac{\bar{x}}{1-\delta_{j}}\left(1+\frac{1}{\bar{z}\gamma_{j}}\right)$
(152)

to satisfy $|\bar{\phi}_{j}|\ll 1$. We will check this latter assumption at
the end of the calculation, but moving forward with it, we find

$\frac{\displaystyle L_{n_{j}-1}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}\left(\bar{\phi}_{j}\right)}=1+\bar{\phi}_{j}+O(\bar{\phi}_{j}^{2})=1-\bar{x}\delta_{j}^{-1}\left(1+\frac{1}{\bar{z}\gamma_{j}}\right)+O(\delta_{j}^{-2}).$
(153)

Inserting Eq.(153) into Eq.(148) and summing over $j$, we obtain

$\displaystyle\langle m\rangle$
$\displaystyle=\sum_{j=1}^{D}\frac{n_{j}}{1-\delta_{j}^{-1}}\frac{\displaystyle
L_{n_{j}-1}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}\left(\bar{\phi}_{j}\right)}$
$\displaystyle=N+\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}-\bar{x}\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}\left(1+\frac{1}{\bar{z}\gamma_{j}}\right)+O(\delta^{-2}).$
(154)

Using Eq.(151), we can write this result as

$\displaystyle\langle m\rangle$
$\displaystyle=N-1+\bar{x}-\bar{x}\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}\left(1+\frac{1}{\bar{z}\gamma_{j}}\right)+O(\delta^{-1}).$
(155)

Therefore, we see that Eq.(155) reproduces Eq.(147) if

$1=\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}\left(1+\frac{1}{\bar{z}\gamma_{j}}\right).$
(156)

With $\langle k\rangle=N-1+O(\delta^{-1})$ and $\bar{z}=N-\langle k\rangle$,
we find $\bar{z}=1+O(\delta^{-1})$, thus giving us the final thermal condition

$1=\sum_{j=1}^{D}{n_{j}}\delta_{j}^{-1}\left(1+\frac{1}{\gamma_{j}}\right)+O(\delta^{-2}).$
(157)

Now we will check the assumption of $|\bar{\phi}_{j}|\ll 1$ for consistency:
From Eq.(157), we can infer that $\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}<1$.
Therefore $\bar{x}$ defined in Eq.(151) is $O(1)$. Also from Eq.(157), we can
infer that, in the large number limit (i.e., $n_{j}\gg 1$), we must have
$\delta_{j}^{-1}(1+1/\gamma_{j})\ll 1$. Thus, we see that
$\bar{\phi}_{j}=\bar{x}(1+{1}/{\bar{z}\gamma_{j}})/(1-\delta_{j})\simeq-\bar{x}\delta_{j}^{-1}(1+1/\gamma_{j})+O(\delta_{j}^{-2})$
can indeed be taken to satisfy $|\bar{\phi}_{j}|\ll 1$.

## References

  * [dM13] Pierre Rémond de Montmort.  Essay d’analyse sur les jeux de hazard… J. Quillau, 1713.
  * [EBC+18] Marie Eubelen, Naguissa Bostaille, Pauline Cabochette, Anne Gauquier, Patricia Tebabi, Andra C Dumitru, Melanie Koehler, Philipp Gut, David Alsteens, Didier YR Stainier, et al.  A molecular mechanism for wnt ligand-specific signaling.  Science, 361(6403):eaat1178, 2018.
  * [EG76] Shimon Even and J Gillis.  Derangements and laguerre polynomials.  Mathematical Proceedings of the Cambridge Philosophical Society, 79:135–143, 1976.
  * [JAVH20] Inga Jarmoskaite, Ishraq AlSadhan, Pavanapuresan P Vaidyanathan, and Daniel Herschlag.  How to measure and evaluate binding affinities.  Elife, 9:e57264, 2020.
  * [JLD10] Iain G Johnston, Ard A Louis, and Jonathan PK Doye.  Modelling the self-assembly of virus capsids.  Journal of Physics: Condensed Matter, 22(10):104101, 2010.
  * [Kra06] Werner Krauth.  Statistical mechanics: algorithms and computations, volume 13.  OUP Oxford, 2006.
  * [MHSB99] Bettina Malnic, Junzo Hirono, Takaaki Sato, and Linda B Buck.  Combinatorial receptor codes for odors.  Cell, 96(5):713–723, 1999.
  * [Mil13] Ron Milo.  What is the total number of protein molecules per cell volume? a call to rethink some published values.  BioEssays, 35(12):1050–1055, 2013.
  * [Nel04] Philip Nelson.  Biological physics: Energy, information life, chapter 8, pages 315–332.  WH Freeman New York, 2004.
  * [PH15] Jason D Perlmutter and Michael F Hagan.  Mechanisms of virus assembly.  Annual review of physical chemistry, 66:217–239, 2015.
  * [PTKG12] Rob Phillips, Julie Theriot, Jane Kondev, and Hernan Garcia.  Physical biology of the cell.  Garland Science, 2012.
  * [RP08] Valerica Raicu and Aurel Popescu.  Integrated molecular and cellular biophysics, volume 10.  Springer, 2008.
  * [RSK04] Bryan L Roth, Douglas J Sheffler, and Wesley K Kroeze.  Magic shotguns versus magic bullets: selectively non-selective drugs for mood disorders and schizophrenia.  Nature reviews Drug discovery, 3(4):353–359, 2004.
  * [SML+20] Christina J Su, Arvind Murugan, James M Linton, Akshay Yeluri, Justin Bois, Heidi Klumpe, Yaron E Antebi, and Michael B Elowitz.  Ligand-receptor promiscuity enables cellular addressing.  bioRxiv, 2020.
  * [Wil18] Mobolaji Williams.  Permutation glass.  Physical Review E, 97(1):012139, 2018.
  * [Wil19] Mobolaji Williams.  Self-assembly of a dimer system.  Physical Review E, 99(4):042133, 2019.