# Derangement model of ligand-receptor binding

Mobolaji Williams  School of Engineering and Applied Sciences,
Harvard University, Cambridge, MA 02138<EMAIL_ADDRESS>
(August 8, 2022)

###### Abstract

We introduce a derangement model of ligand-receptor binding that allows us to
quantitatively frame the question ”How can ligands seek out and bind to their
optimal receptor sites in a sea of other competing ligands and suboptimal
receptor sites?” To answer the question, we first derive a formula to count
the number of partial generalized derangements in a list, thus extending the
derangement result of Gillis and Even. We then compute the general partition
function for the ligand-receptor system and derive the equilibrium expressions
for the average number of bound ligands and the average number of optimally
bound ligands. A visual model of squares assembling onto a grid allows us to
easily identify fully optimal bound states. Equilibrium simulations of the
system reveal its extremes to be one of two types, qualitatively distinguished
by whether optimal ligand-receptor binding is the dominant form of binding at
all temperatures and quantitatively distinguished by the relative values of
two critical temperatures. One of those system types (termed ”search-limited,”
as it was in previous work) does not exhibit kinetic traps and we thus infer
that biomolecular systems where optimal ligand-receptor binding is
functionally important are likely to be search-limited.

Keywords: Derangements, Laguerre Polynomials, Statistical Physics, Ligands and
Receptors, Assembly
MSC codes: 92C05, 82B23, 92-10

###### Contents

  1. 1 Introduction
  2. 2 Partial Generalized Derangements
  3. 3 General Partition Function
    1. 3.1 Gendered Dimer-System Assembly
  4. 4 Large $N$ Limits of Special and General Cases
    1. 4.1 Simple Binding Model
    2. 4.2 Derangement-Only Model
    3. 4.3 General Case
  5. 5 Simulation Comparison
    1. 5.1 Simple Binding Model Simulation $(\delta_{i}=1)$
    2. 5.2 Derangement-Only Model Simulation $(\gamma_{i}\to\infty)$
    3. 5.3 General Case Simulation
  6. 6 Search and Combinatorics Limited Systems
    1. 6.1 Hints at Non-Equilibrium Behavior
  7. 7 Biophysical Implications
  8. 8 Discussion
  9. 9 Conclusion
  10. 10 Acknowledgements
  11. 11 Author’s Statements
  12. 12 Supplementary Code
  13. A Simulations
    1. A.1 Microstate Definition
    2. A.2 Transition Probability
    3. A.3 Transition Types and Examples
  14. B Consistency Checks for $B_{\boldsymbol{r},\boldsymbol{k}}$
    1. B.1 Checking Eq.(8)
    2. B.2 Checking equivalence between Eq.(10) and Eq.(11)
  15. C Deriving General Partition Function
  16. D Laguerre Polynomial Argument Summation Identity Eq.(133)
  17. E Gendered Dimer Assembly Equilibrium Conditions
  18. F Derivation of Eq.(82)

## 1 Introduction

The interaction between membrane receptors and extracellular ligands is the
starting point for many cell-signaling pathways [SML+20]. Given the intricacy
of these pathways, one might think that the initiating ligand-receptor
interaction needs to be “highly specific” (i.e., one ligand type only binds to
one receptor type). But work over the past two decades suggests the opposite:
The specificity of the resulting processes requires such a precise code that
only a combinatorial one, which makes use of various combinations of a finite
number of inputs, can achieve it. This was found in the case of olfactory
receptors [MHSB99] where different receptors recognized different combinations
of ligands. Also, polypharmacology, a recent branch of drug design referring
to creating ligands that act on multiple target receptors, has been found to
be necessary for treating complex diseases such as schizophrenia [RSK04].
Others have found that having many-to-many interactions between ligands and
receptors promotes an increased diversity in range of responses for signaling
pathways [SML+20].

All of these contexts for ligand-receptor interactions allow us to conceive of
a stripped down model of the extracellular medium as one where receptors of
various types and copy numbers exist on a cell surface surrounded by ligands
of various types and copy numbers. Due to the many-to-many interactions
between ligands and receptors (also called, ”multi-specific” or ”promiscuous”
binding), in the most general case such a system exhibits bindings featuring
all combinations of receptors and ligands. Still, to provide a reference point
for the affinities, we can highlight, for each ligand type, a single binding
between receptor and ligand that is the strongest for that ligand. We can term
such interactions as ”optimal” to distinguish them from other interactions.

Such a representation of the ligand-receptor system presents us with a
question: How do the various binding affinities between specific receptors and
ligands affect the global binding properties (e.g., total number of bound
ligands, total number of optimally bound ligands, temperature at which optimal
binding occurs, etc.) of the entire ligand-receptor system?

The binding of ligands to their optimal receptor sites presents both a
combinatorial and a kinetic problem. In a game of musical chairs, a person
sitting in a single chair constrains the chairs available to the remaining
people and thus affects the states the remaining people can occupy. Similarly,
a ligand occupying one receptor site affects the receptor sites available to
other ligands and thus changes the combinatorial count of the possible number
of configurations of those receptors. But for optimal ligand-receptor binding,
ligands not only need to strongly bind to their correct sites; they must also
find these sites. Rather than musical chairs, this aspect of the problem is
more like a game of tag where the size of the environment and the number of
other players determines how easily one person tags another. Similarly, the
size of and the number of receptors in the ligand-receptor system affects how
easily ligands can find receptors.

Thus the combinatorial subproblem for ligand-receptor binding concerns how
ligands can arrange themselves so that each one attaches to its optimal
receptor site, and the kinetic subproblem concerns how ligands can find these
optimal receptor sites in the volume they occupy. The two subproblems together
are also the prototypical definitions of a self-assembly process in which
initially distanced units must come together and combine in the correct
ordered configuration through a thermal system’s unforced evolution towards a
free-energy minimum [Nel04, JLD10, PH15]. There are few analytically tractable
model archetypes that can treat the respective influences of combinatorics and
kinetics on such processes. This work aims to propose such an archetype.

There are some well known approaches to modeling ligand-receptor binding. Most
well known is the law of mass action which has often been applied to ligand-
receptor binding since it provides a coarse-grained framework to model how
affinities affect bound concentrations (see chapter five of [RP08] for a
summary). However, this approach does not take into account the competition
between ligands that can place additional limits on the achievement of
specific types of binding.

A combinatorial model of ligand-receptor binding was presented in chapter six
of [PTKG12]. There the authors considered a collection of identical ligands in
a grid-like space that contained a single receptor. The main combinatorial
task in that analysis was determining the number of ways to arrange the
ligands amongst the spatial grid for both bound and unbound configurations.
From the answer, the authors computed the system partition function and the
ligand-receptor binding concentration. However, the model did not consider
different ligand species existing in the same environment and thus did not
account for the combinatorial competition between various species in real
systems.

In the present work, combinatorics is incorporated through the finiteness of
the number of different particle species in the system and the resulting
finiteness of the possible number of ligand-receptor interactions. Also, by
considering multiple ligand-types each of different copy numbers and with
distinct ”promiscuous” (or multi-specific) binding affinities to a similarly
diverse set of receptor sites, the total possible set of combinatorial
bindings better approaches that of a real system. The net consequence of these
assumptions is to introduce combinatorial competition into the equilibrium
statistical physics that define the system, making finite number effects
particularly important in describing thermal properties.

The general system we consider is shown in Fig. 1. Say we have a system of
ligands and receptors existing in the extracellular medium. The ligands come
in multiple copies as do the receptors, and, as is consistent with the multi-
specificity of some real ligand systems, we assume all ligands have the
ability to bind to all receptors. However, we will also assume that each
ligand type has a specific optimal binding with a particular receptor type.
This latter assumption will provide us with an additional order parameter with
which we can define our system.

Figure 1: Example ligand-receptor system. The figure displays three different
types of ligands and three corresponding types of receptors. All ligands can
bind to all receptors, but each receptor has an optimal binding with a
specific receptor type. The ligands not bound to any receptor are free and
exist in the volume of the system. In this work we derive the conditions under
which all ligands can bind to their optimal receptor sites.

There are some basic questions that a corresponding model for this system
should be able to answer: How does the average number of bound ligands of
various types depend on the system’s binding affinities? How does the average
number of optimally bound ligands of various types depend on the system’s
binding affinities? What thermal conditions define a system in which all
ligands are optimally bound? For what properties of the affinities are such
conditions even feasible? Are there ways we can categorize these systems so as
to determine a priori from affinity properties what the expected binding
behavior should be?

We answer these questions in the subsequent sections. In Section 2, we
introduce the main combinatorial problem that underlies the binding of
multiple ligand species to multiple receptor species. In Section 3, we use the
solution to this problem to compute the partition function for the system and
show that the result generalizes a special case derived in [Wil19]. In Section
4, we consider the large particle-number limit of the partition function for
two limiting cases and the general case. The limiting cases help us build the
intuition relevant to understanding the conditions that define fully optimal
binding in the general case. In each case, we derive expressions for the
average number of bound ligands of each species and the average number of
optimally bound ligands, as long as the relevant quantity is not trivially
constrained by the case itself. In Section 5, we introduce an image on a grid
to give a visual handle on the various limiting cases of the system, and we
simulate the grid images to affirm that the analytical results accurately
predict the average grid state at various temperatures. In Section 6, we note
that the temperature curves for the average number of bound ligands and the
average number of optimally bound ligands have distinct limiting behaviors
contingent on how model parameters vary with one another. We explore these
distinct limiting behaviors through simulations and argue that one limiting
behavior is associated with kinetic traps. In Section 7, we return to the
system that motivated our analyses and discuss the biophysical implications of
the results. In Section 8, we conclude by considering ways to extend the
general model.

## 2 Partial Generalized Derangements

Our ultimate goal is to model the equilibrium thermodynamics of systems of the
kind shown in Fig. 1. Achieving this amounts to computing a partition function
and then using the partition function to find order parameters, but first we
need to solve the combinatorial problem at the heart of this system.

We recall that the derangement of a list is a rearrangement of that list such
no element is in its original position. The formula for the number of
derangements of a list with $N$ unique elements was first obtained by Pierre
Mortmort and Nicholas Bernoulli in the the early 18th century [dM13]:

$d_{N}=\sum_{j=0}^{N}\binom{N}{j}(-1)^{j}(N-j)!.$ (1)

More than 200 years later, Gillis and Even derived the generalization to this
result for the case where elements occur with multiple copy number [EG76].
They showed that the number of ways to completely derange an ordered list with
$n_{1}$ elements of type 1, $n_{2}$ elements of type $2$, $\ldots$, and
$n_{D}$ elements of type $D$ (where $D\leq N$) is

$G_{\boldsymbol{n}}=\int^{\infty}_{0}dx\,e^{-x}\prod_{k=1}^{D}(-1)^{n_{k}}L_{n_{k}}(x),$
(2)

where $\boldsymbol{n}\equiv(n_{1},n_{2},\ldots,n_{D})$ and $L_{n}(x)$ is the
Laguerre polynomial defined as

$L_{n}(x)=\sum_{j=0}^{n}\binom{n}{k}\frac{(-1)^{j}}{j!}x^{j}.$ (3)

We will call Gillis and Even’s result the ”generalized derangement result.”
For this work, we want to obtain a further generalization to the generalized
derangement result to the case where not necessarily all elements of an
initial list are included in a rearrangement. Finding this generalization
would allow us to model a system in which ligands can exist both on and off
receptor sites.

The primary problem we need to solve is as follows:

> We have $r_{1}$ elements of type $1$, $r_{2}$ elements of type $2$,
> $\ldots$, and $r_{D}$ elements of type $D$, all of which are arranged in an
> initial list. All elements are then removed from the list. What is the
> number of ways that we can choose and arrange $k_{1}\leq r_{1}$ elements of
> type $1$, $k_{2}\leq r_{2}$ elements of type $r_{2}$, $\ldots$, and
> $k_{D}\leq r_{D}$ elements of type $D$ such that none of the elements has
> the same position as it has in the original list?

We call the answer to this question the ”partial generalized derangement
result,” given that we are considering derangements of partial collections of
the total set of elements with repeats. The resulting quantity will be denoted
$B_{\boldsymbol{r},\boldsymbol{k}}$ where
$\boldsymbol{r}=(r_{1},r_{2},\ldots,r_{D})$ and
$\boldsymbol{k}=(k_{1},k_{2},\ldots,k_{D})$, and we will obtain an explicit
expression for it by reasoning according to the principal of inclusion and
exclusion.

To apply the principal of inclusion and exclusion in the desired case, it is
helpful to first review it in the simpler case of Eq.(1). With the summation
index $j$ denoting the number of elements that are fixed in their original
positions, the factor $\binom{N}{j}$ is the number of ways to choose $j$ fixed
elements out of $N$ possible elements. The factor $(-1)^{j}$ is the common
principle of inclusion and exclusion factor that leads sets of ”correct
position” elements to be alternately subtracted from and added to the first
term of $N!$ which is a count of all permutations. The factor $(N-j)!$ counts
the number of ways to arrange the remaining elements given that $j$ are fixed
in their original positions. The end result after summing over all $j$ is a
count of only permutations that do not include any elements in their original
positions.

Thus, there are three essential factors in the summand of Eq.(1): The counting
of the number of ways to arrange elements in their original position; the
principle of inclusion and exclusion factor $(-1)$ for each such original-
position element; and the factor that counts the number of ways to arrange the
remaining elements.

We can define analogous factors for $B_{\boldsymbol{r},\boldsymbol{k}}$ and
use them to write a summation expression for the quantity. The result is

$\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})!}.$
(4)

To understand Eq.(4) we consider how each factor in the summand contributes to
the final expression. The factor $\binom{r_{i}}{j_{i}}$, for $i=1,\ldots,D$,
is the number of ways to fill $j_{i}$ out of the $r_{i}$ positions of type $i$
with their original elements. The factor $(-1)^{j_{1}+\cdots+j_{D}}$ is the
net principal of inclusion and exclusion factor for the $j_{i}$ elements of
type $i$ (for $i$ running from $1$ to $D$) that are in their original
positions. After fixing these positions with their original elements, there
are now $r_{1}-j_{1}+\cdots+r_{D}-j_{D}$ possible positions which we must fill
with $k_{1}-j_{1}+\cdots+k_{D}-j_{D}$ elements. The number of ways to choose
which of these remaining positions to fill is represented by a binomial
factor. The last factor
${(k_{1}-j_{1}+\cdots+k_{D}-j_{D})!}/{(k_{1}-j_{1})!\cdots(k_{D}-j_{D})!}$ is
the number of ways to permute the $k_{1}-j_{1}+\cdots+k_{D}-j_{D}$ elements
amongst the chosen positions divided by factors to correct for the fact that
elements of the same type are identical.

To affirm correctness, we can perform some sanity checks on Eq.(4) to show
that this result is consistent with related ones.

For what follows, it will be most useful to express Eq.(4) as an integral
expression. To do so, we introduce the generalized Laguerre polynomial:

$L_{n}^{(\alpha)}(x)=\sum_{j=0}^{n}\binom{n+\alpha}{n-j}\frac{(-1)^{j}}{j!}x^{j}.$
(5)

Using Eq.(5) and the definition of the Gamma function, we find that Eq.(4) can
be written as

$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),$
(6)

where we defined $\alpha_{i}\equiv n_{i}-k_{i}$.

Figure 2: Microstate of partial generalized derangements. In a partial
generalized derangement, elements occur in multiple copies and partially
occupy deranged positions in a list. The figure shows a microstate for the
system with $D=3$, $r_{1}=6$, $r_{2}=3$, and $r_{3}=5$ with squares,
triangles, and circles associated with $1$, $2$, and $3$ respectively. There
are $k_{1}=3$, $k_{2}=2$, and $k_{3}=3$ elements of the various types in
contact with the lattice. Some of the elements on the lattice are in deranged
positions and some are in correct positions. Taking
$\boldsymbol{m}=(m_{1},m_{2},m_{3})$ to define the vector counting the number
of elements of each type in correct positions, we have $m_{1}=2$, $m_{2}=1$,
and $m_{3}=1$. Given these correct positions, the microstate in this figure
contributes to the count for
$B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}$.

For the first sanity check, we expect that Eq.(6) should reduce to Eq.(2) when
we take $k_{j}=n_{j}$ for all $j$. Namely when we are considering the full
(rather than a partial) set of elements, the partial generalized derangement
result should reduce to the generalized derangement result. Imposing this
equality on Eq.(6) and noting that $L_{j}\equiv L^{(0)}_{j}$, we indeed find
that Eq.(2) is reproduced.

One can show (as was done in the appendix of [Wil19]) that if we have $D$
different elements each of which is associated with a particular site out of
$D$ lattice sites, then the number of ways to select $K\leq D$ elements to
arrange amongst the lattice sites such that none is in its associated site is

$b_{D,K}=\sum_{J=0}^{D}(-1)^{J}\binom{D}{J}\binom{D-J}{K-J}^{2}(K-J)!.$ (7)

Thus, we should be able to show that Eq.(6) reduces to Eq.(7) under the right
conditions. In particular if we take
$\boldsymbol{r}=(1,1,\ldots,1)\equiv\boldsymbol{r}_{0}$ (i.e., we have $D$
unique elements, each of a single copy-number), then the vector
$\boldsymbol{k}$ in Eq.(6) can only have elements of $1$ or $0$, and thus
$\boldsymbol{k}$ defines a particular subset of the total set of elements.
$B_{\boldsymbol{r}_{0},\boldsymbol{k}}$ then represents the number of ways to
completely derange a particular collection of unique elements where the
collection is defined by the vector $\boldsymbol{k}$. In order to find the
total number of ways to completely derange $K$ total elements (i.e., what is
represented in Eq.(7)), we need to sum $B_{\boldsymbol{r}_{0},\boldsymbol{k}}$
over all possible values of $\boldsymbol{k}$ such that $\sum_{j}k_{j}=K$. Thus
the consistency check we must make is

$\sum_{k_{1}=0}^{1}\ldots\sum_{k_{D}=0}^{1}B_{\boldsymbol{r}_{0},\boldsymbol{k}}\delta(K,k_{1}+\cdots+k_{D})=b_{D,K}.$
(8)

It takes more work to demonstrate Eq.(8) (see Appendix B.1), but doing so
affirms that Eq.(6) is consistent with its simpler manifestations.

As a final consistency check, we note that there should be a summation
condition for the total number of ways to select $k_{1}$ elements of type $1$,
$k_{2}$ elements of type $2$, $\ldots$, $k_{D}$ elements of type $D$ to be
placed amongst the available $r_{1}+\cdots+r_{D}$ positions without regard to
whether the elements are placed in their original positions. For the case of
simple derangements Eq.(1), this summation condition is

$N!=\sum_{m=0}^{N}\binom{N}{m}d_{N-m}.$ (9)

Eq.(9) represents the fact that the total number of ways to order $N$ unique
elements is also the number of ways to select $m$ fixed elements and derange
the rest summed over all possible values of $m$. It is straightforward to
check that Eq.(1) satisfies Eq.(9).

Towards finding an analogous summation condition for
$B_{\boldsymbol{r},\boldsymbol{k}}$, we note that
$\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_{i}$ out of $r_{i}$ positions (for
$i=1,\ldots,D$) to contain their original elements while the remaining
$k_{i}-m_{i}$ elements are completely deranged with respect to the
$r_{i}-m_{i}$ remaining original positions of type $i$. If we sum this
quantity over all possible values of $m_{i}$, as in

$I_{\boldsymbol{r},\boldsymbol{k}}\equiv\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}},$
(10)

we should obtain the number of ways to arrange (and not necessarily derange)
$k_{i}\leq r_{i}$ elements of type $i$ for $i=1,\ldots,D$ across a total of
$r_{1}+\ldots+r_{D}$ lattice sites.

Calculating this quantity another way, we note that (including filled and
empty sites) we are technically trying to order a total of
$r_{1}+\cdots+r_{D}$ sites: There are $k_{1}+\cdots+k_{D}$ filled sites and
$r_{1}-k_{1}+\cdots+r_{D}-k_{D}$ empty sites. Consequently there are
$(r_{1}+\cdots+r_{D})!$ ways to order the total collection. Given that the
filled-site elements occur in multiple copies, we must correct for equivalent
orderings by dividing this count by $k_{j}!$ for each element type. Also,
since the empty sites act as an extra ”type” of element, we must also divide
the count by $(r_{1}-k_{1}+\cdots+r_{D}-k_{D})!$, the number of ways to
reorder these empty sites. Thus we should find

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

In Appendix B.2, we show that Eq.(10) produces Eq.(11).

With our combinatorial expression found and consistency affirmed, we can now
work towards building the partition function for the system.

## 3 General Partition Function

We recall that our objective is to study the equilibrium thermodynamics of the
physical system depicted in Fig. 1. The system is one where a fixed set of
ligands can exist as bound or unbound to a collection of receptors. When a
ligand is bound to a receptor, it can be bound either to an optimal receptor
or to a suboptimal receptor. To study the thermodynamics of such a system, we
needed to compute a combinatorial factor that counts the number of ways
ligands can be bound to receptor sites where some of these bindings are
suboptimal. Having computed this quantity in Sec. 2, we can now use what we
found to calculate the partition function.

However, before using this derangement formalism, we will begin with minimal
assumptions and write the most general expression possible for the partition
function of the system. The intent in starting here is to show the
intractability of the most general form of the partition function and thereby
demonstrate the analytical benefits afforded by considering derangements from
a pre-defined sequence.

We start by defining numerical quantities in the system. Say that we have $D$
different types of ligands and a corresponding set of $D$ different types of
receptors. A ligand type and a receptor type are labeled with $i$ for
$i=1,\ldots,D$. The ligand of type $i$ has $n_{i}$ copies in the system, and
the receptor of type $i$ has $r_{i}$ copies in the system. Each ligand can
either be bound to a receptor or be unbound and free to move in the space
surrounding the receptor sites. There are $N_{R}\equiv r_{1}+\cdots+r_{D}$
total receptors and each of the $N_{L}\equiv n_{1}+\cdots+n_{D}$ ligands can
bind to any one of them, provided there is an available binding site.

In the most general theoretical formulation of the problem, we can represent
the system microstate by a matrix ${\cal C}$ with elements ${\cal C}_{i,j}$
defined as

${\cal C}_{i,j}=\text{$\\#$ of bindings between ligand of type $i$ and
receptor of type $j$}.$ (12)

If we can specify all elements of the $D\times D$ matrix ${\cal C}$ then we
have completely defined the system. Given our counts for the total number of
ligands and receptors of each type, there are only three constraints on the
elements ${\cal C}_{i,j}$: Each element must be an integer,
$\sum_{j=1}^{D}{\cal C}_{i,j}\leq n_{i}$, and $\sum_{i=1}^{D}{\cal
C}_{i,j}\leq r_{j}$.

Next, we ask how we can incorporate binding parameters to represent the way
energy affects the probability of a microstate. We will take $Q_{i,j}^{B}$ to
be the single-particle partition function for a ligand of type $i$ that is
bound to a receptor of type $j$. Therefore, the multi-particle partition
function for all ligands of type $i$ that are bound to receptors of type $j$
is $(Q^{B}_{i,j})^{{\cal C}_{i,j}}$. There is no factorial correction in this
expression because ligands that are bound to specific receptor sites are
distinguishable by virtue of the distinguishability of the receptor sites
themselves. Conversely, taking $Q_{i}^{F}$ to be the single-particle partition
function for an unbound ligand of type $i$, and given that
$n_{i}-\sum_{j=1}^{D}{\cal C}_{i,j}$ is the total number of unbound ligands of
type $i$, the multi-particle partition function for the ligand of type $i$
that is unbound is ${\left(Q^{F}_{i}\right)^{n_{i}-\sum_{j=1}^{D}{\cal
C}_{i,j}}}/{(n_{i}-\sum_{j=1}^{D}{\cal C}_{i,j})!}$ where the factorial
correction is because these ligands are in free space. Putting the pieces
together, and including appropriate products to account for various ligand and
receptor types, we find that the general partition function for this system is

${\cal Z}_{\boldsymbol{n},\boldsymbol{r}}=\sum_{\\{{\cal
C}_{i,j}\\}}\prod_{i=1}^{D}\frac{\left(Q^{F}_{i}\right)^{n_{i}-\sum_{j=1}^{D}{\cal
C}_{i,j}}}{(n_{i}-\sum_{j=1}^{D}{\cal
C}_{i,j})!}\,\prod_{j=1}^{D}(Q^{B}_{i,j})^{{\cal C}_{i,j}},$ (13)

where $\boldsymbol{n}\equiv(n_{1},\ldots,n_{D})$ and
$\boldsymbol{r}=(r_{1},\ldots,r_{D})$. The product over $j$ is the product
over all types of receptors to which ligand $i$ is bound. The product over $i$
is the product over all types of ligands, with the associated factors
representing free and bound ligands. The external summation is a sum over all
possible values of ${\cal C}_{i,j}$ according to the constraints of the
physical system (i.e., ${\cal C}_{i,j}$ is an integer, $\sum_{j=1}^{D}{\cal
C}_{i,j}\leq n_{i}$ and $\sum_{i=1}^{D}{\cal C}_{i,j}\leq r_{j}$).

With Eq.(13), we have our exact partition function for the system. However, it
is unclear how to make this exact expression analytically useful. Ostensibly
we need to enumerate and then sum over all possible matrices ${\cal C}$, but
this is unfeasible given the number of microstates associated with even
simpler contact matrices. For the products, simplifications often occur when
products can be turned into summations, but the matrix nature of the factors
in Eq.(13) seems to preclude this conversion. Therefore, we cannot compute
analytical expressions for observables from Eq.(13) since such quantities
depend on tractable calculations of the partition function. Because of these
reasons, we need to make some simplifying assumptions to make progress.

In our system, each ligand of type $i$ will bind to a receptor of type $j$
with a binding Boltzmann factor $Q_{i,j}^{B}$. The matrix nature of this
expression seems to be the principal complication in our partition function
since it prevents us from converting the product over $j$ into a sum in the
power. Therefore, to simplify this expression we will make an assumption that
reduces the dimensionality of the parameters space of this matrix.

While still assuming that any ligand can bind to any receptor, we will also
assume that every ligand of type $i$ has the same binding affinity to every
type of receptor except to a receptor of type $i$ to which the ligand of type
$i$ binds with an additional energy of binding of $\Delta_{i}$ (with
$\Delta_{i}\geq 0$). Mathematically, we can encode this assumption into the
model by making the transformation

$Q_{i,j}^{B}\quad\longrightarrow\quad(Q^{B}_{i})e^{\beta\Delta_{i}\delta_{i,j}},$
(14)

where $\delta_{i,j}$ is the Kronecker delta, and $\beta=1/k_{B}T$ for a system
at temperature $T$. Essentially, Eq.(14) asserts that all ligands of type $i$
that are bound to a receptor have the same single-particle partition function
$Q^{B}_{i}$ except for the type $i$ ligands that are bound to type $i$
receptors which acquire an additional Boltzmann factor $e^{\beta\Delta_{i}}$.
We call the latter such bindings ”optimal” or ”correct” bindings; all other
bindings are termed ”sub-optimal.”

One benefit of Eq.(14) is that it reduces the number of parameters that are
needed to define the system: Rather than have $N^{2}$ binding parameters
defined by the elements of $Q_{i,j}^{B}$, we have $2N$ parameters from
$(Q^{B}_{i})$ and $\Delta_{i}$ together, a reduction that makes our modeling
more tractable.

Regarding the physical motivation of this assumption, Eq.(14) closely
approximates the ligand-receptor interactions of signaling pathways such as
the Wnt-Fz pathway which exhibits specific binding for the initiation of
signals (i.e., a nonzero $\Delta_{i}$) while still having ligands with
promiscuous interactions with receptors [EBC+18].

How does Eq.(14) change our partition function? Starting from Eq.(13) and
incorporating the transformation Eq.(14), we obtain

$\displaystyle{\cal Z}_{\boldsymbol{n},\boldsymbol{r}}$
$\displaystyle=\sum_{\\{{\cal
C}_{i,j}\\}}\prod_{i=1}^{D}\frac{\left(Q^{F}_{i}\right)^{n_{i}-\sum_{j=1}^{D}{\cal
C}_{i,j}}}{(n_{i}-\sum_{j=1}^{D}{\cal
C}_{i,j})!}\,(Q^{B}_{i})^{\sum_{j=1}^{D}{\cal
C}_{i,j}}e^{\beta\Delta_{i}\sum_{j=1}^{D}{\cal C}_{i,j}\delta_{i,j}}$ (15)

From, here we make a change of variables motivated by the quantities that
appear in the expression. We define $k_{i}=\sum_{j=1}^{D}{\cal C}_{i,j}$ as
the total number of bound (optimally or not) ligands of type $i$ and
$m_{i}=\sum_{j=1}^{D}{\cal C}_{i,j}\delta_{j,i}$ as the total number of
optimally bound ligands of type $i$. Eq.(15) then becomes

$\displaystyle{\cal Z}_{\boldsymbol{n},\boldsymbol{r}}$
$\displaystyle=\sum_{\boldsymbol{k}}\sum_{\boldsymbol{m}}\Omega(\boldsymbol{k},\boldsymbol{r},\boldsymbol{m})\prod_{i=1}^{D}\frac{\left(Q^{F}_{i}\right)^{n_{i}-k_{i}}}{(n_{i}-k_{i})!}\,(Q^{B}_{i})^{k_{i}}e^{\beta\Delta_{i}m_{i}}$
(16)

where we defined the summations as

$\sum_{\boldsymbol{k}}\equiv\prod_{i=1}^{D}\sum_{k_{i}=0}^{n_{i}},\qquad\sum_{\boldsymbol{m}}\equiv\prod_{i=1}^{D}\sum_{m_{i}=0}^{k_{i}},$
(17)

and the combinatorial factor as

$\Omega(\boldsymbol{k},\boldsymbol{r},\boldsymbol{m})=\sum_{\\{{\cal
C}_{i,j}\\}}\prod_{i=1}^{D}\delta\left(k_{i},\textstyle\sum_{j=1}^{D}{\cal
C}_{i,j}\right)\delta\left(m_{i},\textstyle\sum_{j=1}^{D}{\cal
C}_{i,j}\delta_{i,j}\right)$ (18)

resulting from the change of summation variables. Qualitatively
$\Omega(\boldsymbol{k},\boldsymbol{r},\boldsymbol{m})$ is the number of ways
to fill $m_{i}$ out of $r_{i}$ receptor sites (for $i=1,\ldots,D$) with their
optimal binding partner ligands while having the remaining $k_{i}-m_{i}$ bound
ligands (again for $i=1,\ldots,D$) not in their associated optimal site.
Eq.(15) is already an improvement over Eq.(13) since there are no longer any
matrices as factors and our summations are over specific integer-valued
variables (i.e., $k_{i}$ and $m_{i}$) rather than elements of a set (i.e.,
$\\{{\cal C}_{i,j}\\}$). However, it seems that in order to compute Eq.(16),
we would need to compute
$\Omega(\boldsymbol{k},\boldsymbol{r},\boldsymbol{m})$ which according to
Eq.(17) seems to require a summation over the elements of the aforementioned
set. Fortunately this is not the case: We already have the necessary
quantities to compute this factor. Given the qualitative definition of
$\Omega$, we can assert that

$\Omega(\boldsymbol{k},\boldsymbol{r},\boldsymbol{m})=B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}\prod_{i=1}^{D}\binom{r_{i}}{m_{i}},$
(19)

where $B_{\boldsymbol{n},\boldsymbol{k}}$ is defined in Eq.(6). On the right-
hand side of Eq.(19), the factor
$B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}$ represents
the number of ways to select and arrange $k_{i}-m_{i}$ ligands of type $i$
(for $i=1,\ldots,D$) across a total set of $r_{1}-m_{1}+\ldots r_{D}-m_{D}$
receptors such that no ligand of type $i$ is bound to one of its $r_{i}-m_{i}$
optimal receptors. The factors
$\binom{r_{1}}{m_{1}}\cdots\binom{n_{D}}{m_{D}}$ count the number of ways to
choose $m_{i}$ receptors from the $r_{i}$ possible receptors for
$i=1,\ldots,D$ to be occupied by optimal-binding partner ligands. Thus, we see
that the product of
$B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}$ and
$\binom{r_{1}}{m_{1}}\cdots\binom{n_{D}}{m_{D}}$ is indeed the number of ways
to fill $m_{i}$ out of $r_{i}$ receptor sites (for $i=1,\ldots,D$) with their
optimal binding partner ligands while having the remaining $k_{i}-m_{i}$ bound
ligands (again for $i=1,\ldots,D$) not in their associated optimal sites.

With the observation that $\Omega$ can be computed from Eq.(19), we can now
obtain a simpler expression for Eq.(16): We have

$\displaystyle{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}=\sum_{\boldsymbol{k}}\sum_{\boldsymbol{m}}B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{k}-\boldsymbol{m}}\prod_{i=1}^{D}\binom{r_{i}}{m_{i}}e^{\beta
m_{i}\Delta_{i}}(Q_{i}^{B})^{k_{i}}\frac{(Q_{i}^{F})^{n_{i}-k_{i}}}{(n_{i}-k_{i})!},$
(20)

We recall that in the summands of Eq.(20), $k_{i}\leq n_{i}$ represents the
total number of ligands of type $i$ that are bound to receptors, and
$m_{i}\leq k_{i}$ represents the total number of ligands of type $i$ that are
optimally bound to receptors.

For notational simplicity, we will define some additional constants. We define

$c_{\boldsymbol{n}}\equiv\prod_{i=1}^{D}{(Q_{i}^{F})^{n_{i}}},\qquad\delta_{i}\equiv
e^{\beta\Delta_{i}},\qquad\gamma_{i}\equiv\frac{Q_{i}^{B}}{Q_{i}^{F}},$ (21)

thus giving us

$\displaystyle{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}=c_{\boldsymbol{n}}\sum_{\boldsymbol{k}}\sum_{\boldsymbol{m}}B_{\boldsymbol{r}-\boldsymbol{m},\,\boldsymbol{k}-\boldsymbol{m}}\prod_{i=1}^{D}\binom{r_{i}}{m_{i}}\frac{1}{(n_{i}-k_{i})!}\,\delta_{i}^{m_{i}}\gamma_{i}^{k_{i}}.$
(22)

Given that pre-factors do not affect physical predictions in canonical
partition functions, Eq.(22) reveals that it is only the ratios of our single-
particle partition functions that are thermodynamically relevant. This result
makes sense given that only free-energy differences (i.e., logarithms of
partition function ratios) should affect the physics of a system. Thus,
without loss of generality, we can impose $Q_{i}^{F}=1$ for all $i$ under the
assumption that the thermal dependence of each $Q_{i}^{F}$ can be absorbed
into a redefinition of $\gamma_{i}$ and $\delta_{i}$ with no change in the
physical implications of Eq.(22). With this imposition we have

$c_{\boldsymbol{n}}=1\qquad\text{[free-particle partition function
normalization]}$ (23)

Moving forward, we recognize that the partition function becomes more
analytically useful to us if we can replace the discrete summation with an
integral since integrals, unlike discrete summations, are more amenable to the
methods of analysis. To do so we make use of the integral form of
$B_{\boldsymbol{n},\boldsymbol{k}}$ in Eq.(6) and a few Laguerre polynomial
identities. After some work (see Appendix C), we obtain

${\cal Z}_{\boldsymbol{n},\boldsymbol{r}}=\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).$
(24)

We can write this result in a more mathematically useful form by expressing
the integrand as the exponential of a potential function:

$\displaystyle{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta},\boldsymbol{\gamma})=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z}\int^{\infty}_{0}dx\,\exp\big{[}{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\big{]},$
(25)

where $\Gamma$ is a closed contour about the origin in the complex plane and

$\displaystyle{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})$
$\displaystyle=z-x+(N_{R}-N_{L})\ln(x/z)$
$\displaystyle\qquad\qquad+\sum_{j=1}^{D}\ln\left[(\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)\right],$
(26)

with $L^{(\alpha)}_{n}(x)$ the $n$th generalized Laguerre polynomial, and
$N_{R}\equiv\sum_{j=1}^{D}r_{j}$ and $N_{L}\equiv\sum_{j=1}^{D}n_{j}$,
respectively, the total number of receptors and total number of ligands in the
system. Eq.(25) provides the starting point for our thermal equilibrium
analysis. But first we will derive expressions for the order parameters
written in terms of this partition function.

From Eq.(22), we can derive expressions for the two main observables of the
system. The average number of bound ligands and the average number of
optimally-bound ligands are, respectively,

$\langle k\rangle=\sum_{j=1}^{D}\langle
k_{j}\rangle=\sum_{j=1}^{D}\gamma_{j}\frac{\partial}{\partial\gamma_{j}}\ln{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}},\qquad\langle
m\rangle=\sum_{j=1}^{D}\langle
m_{j}\rangle=\sum_{j=1}^{D}\delta_{j}\frac{\partial}{\partial\delta_{j}}\ln{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}.$ (27)

We can use the second equation in Eq.(27) to write an alternative expression
for $\langle m\rangle$. For the function
$f_{n}^{(\alpha)}(x;q)=(q-1)^{n}L_{n}^{(\alpha)}(x/(1-q))$, the identity
$u\,\partial_{u}L_{n}^{(\alpha)}(u)=nL_{n}^{(\alpha)}(u)-(n+\alpha)L_{n-1}^{(\alpha)}(u)$
allows us to prove

$\partial_{q}f_{n}^{(\alpha)}(x;q)=(n+\alpha)f_{n-1}^{(\alpha)}(x;q).$ (28)

We can then show

$\displaystyle\frac{\partial}{\partial\delta_{i}}\left[(\gamma_{i}(\delta_{i}-1))^{n_{i}}L_{n_{i}}^{(r_{i}-n_{i})}\left(\frac{x(z\gamma_{i}+1)}{z\gamma_{i}(1-\delta_{i})}\right)\right]$
$\displaystyle=r_{i}\gamma_{i}(\gamma_{i}(\delta_{i}-1))^{n_{i}-1}L_{n_{i}-1}^{(r_{i}-n_{i})}\left(\frac{x(z\gamma_{i}+1)}{z\gamma_{i}(1-\delta_{i})}\right).$
(29)

Thus from Eq.(27), we have

$\langle m\rangle=\sum_{j=1}^{D}r_{j}\gamma_{j}\delta_{j}\frac{{\cal
Z}_{\boldsymbol{n}_{j},\boldsymbol{r}_{j}}(\boldsymbol{\delta},\boldsymbol{\gamma})}{{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta},\boldsymbol{\gamma})},$
(30)

where $\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})$. The
vector $\boldsymbol{r}_{j}$ is defined similarly. Although the term $r_{i}-1$
does not appear in Eq.(29), we had to introduce $\boldsymbol{r}_{j}$ into the
expression Eq.(30) to ensure that $r_{j}-n_{j}$ in Eq.(24) remained unchanged
when we replaced $\boldsymbol{n}$ with $\boldsymbol{n}_{j}$. Eq.(30) could
also have been derived from Eq.(22) by differentiating with respect to
$\delta_{j}$ and making the dummy variable replacement
$k_{j}^{\prime}=k_{j}-1$ and $m^{\prime}_{j}=m_{j}-1$. There is no simplified
expression for $\langle k\rangle$ analogous to Eq.(30).

As is common for partition functions written as integrals, approximating the
partition function by the maximum (or, in the case of complex values,
stationary) value of its integrand allows us to derive more tractable
expressions for the equilibrium conditions. Before we pursue these conditions,
we will show how Eq.(25) is a generalization of a result established in a
previous paper. The purpose of establishing such a generalization is to
extrapolate some of the physical results explored in that previous paper to
this more complex case.

###  3.1 Gendered Dimer-System Assembly

Figure 3: System for gendered dimer assembly: The general system studied in
this work for the case of $n_{i}=1$. For this case, there is one copy of each
particle species and each particle has a one-to-one correspondence with an
optimal binding site. Also, the binding affinities and optimal-binding
affinity advantages for each particle vary with the particle species.
Consequently, this case is a slight generalization of the case of dimer
assembly with fixed binding sites treated in [Wil19].

In the appendix of [Wil19], we analyzed a physical system termed ”gendered
dimer assembly.” We recall that this refers to a system where there are two
types of particles and where a particle of one type can only bind to a
particle of the other type. By fixing the positions of all particles of one
type, we were able to apply the results to the case of particles binding to a
lattice of possible sites (e.g., ligands binding to receptors). We assumed
each particle type had a single copy and that all particles had the same
binding affinities to the lattice and the same optimal-binding affinities to
their correct sites.

In this section we use the general expression Eq.(25) to derive a slight
generalization of these past results. In [Wil19] we assumed a global binding
affinity and optimal-binding affinity for all types of particles, but here we
will assume that particles’ binding affinities and optimal-binding affinities
vary according to particle type. In effect, in Eq.(25) we will take
$r_{i}=n_{i}=1$ for all $i$ but retain the index dependence of $\delta_{i}$
and $\gamma_{i}$. Since gendered dimer assembly (with one of the ”genders”
fixed in space) is a more specific case of the ligand-receptor binding system
considered in this work, we should find that the equilibrium equations derived
for this special case of Eq.(25) match those found for the gendered dimer
system in [Wil19].

First, imposing the condition $r_{i}=n_{i}=1$ for all $i$ on Eq.(25), we find
the partition function

$Z_{D}(\boldsymbol{\delta},\boldsymbol{\gamma})=\frac{1}{2\pi
i}\oint\frac{dz}{z}\int^{\infty}_{0}dx\,\exp\left[F_{D}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\right],$
(31)

where

$F_{D}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})=z-x+\sum_{i=1}^{D}\ln\Big{(}\gamma_{j}(\delta_{j}-1+x)+\frac{x}{z}\Big{)}.$
(32)

The corresponding average number of bound particles and average number of
optimally bound particles for a particle of type $j$ are the same as what is
given in Eq.(27):

$\langle k_{j}\rangle=\gamma_{j}\frac{\partial}{\partial\gamma_{j}}\ln
Z_{D},\qquad\langle
m_{j}\rangle=\delta_{j}\frac{\partial}{\partial\delta_{j}}\ln Z_{D}.$ (33)

Applying the large $N$ integral approximation (specifically $D\gg 1$ in this
case) to Eq.(31) yields

$Z_{D}(\boldsymbol{\delta},\boldsymbol{\gamma})\simeq\frac{1}{(\bar{z}^{2}\det
H)^{1/2}}\exp\left[F_{D}(\bar{z},\bar{x};\boldsymbol{\delta},\boldsymbol{\gamma})\right],$
(34)

where $H$ is the hessian matrix with second-order derivative components
$H_{\alpha,\beta}=\partial_{\alpha,\beta}F_{D}\big{|}_{z=\bar{z},x=\bar{x}}$
($\alpha,\beta\in\\{x,z\\}$), and $\bar{z}$ and $\bar{x}$ are defined by the
conditions

$0=\partial_{x}F_{D}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\Big{|}_{z=\bar{z}x=\bar{x}},\qquad
0=\partial_{z}F_{D}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\Big{|}_{z=\bar{z},x=\bar{x}}.$
(35)

In order to compute equilibrium conditions for $\langle k_{j}\rangle$ and
$\langle m_{j}\rangle$ from Eq.(34) we first need the conditions for $\bar{x}$
and $\bar{z}$222When applying the saddle-point approximation or Laplace’s
method to an integral, one should check that the second derivative matrix has
the correct stability properties to ensure the validity of the approximation.
In the most general case in this work, this check requires us to compute the
spectral properties of a $2\times 2$ complex matrix. The difficulty in proving
stability in this general scenario has led us to instead use simulations to
heuristically vet the validity of the approximation. Using Eq.(32) and Eq.(35)
to find these conditions, we have, from $\partial_{x}F_{D}=0$ and
$\partial_{z}F_{D}=0$, respectively,

$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})}.$
(36)

Next, computing $\langle k_{j}\rangle$ and $\langle m_{j}\rangle$, we have

$\displaystyle\langle k_{j}\rangle$
$\displaystyle=\gamma_{j}\frac{\partial}{\partial\gamma_{j}}F_{D}=\frac{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})}{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})+\bar{x}}$
(37) $\displaystyle\langle m_{j}\rangle$
$\displaystyle=\delta_{j}\frac{\partial}{\partial\delta_{j}}F_{D}=\frac{\bar{z}\gamma_{j}\delta_{j}}{\bar{z}\gamma_{j}(\delta_{j}-1+\bar{x})+\bar{x}}$
(38)

where, in applying Eq.(33) to Eq.(34), we neglected the exponential pre-factor
in the latter since it is subleading in the $D\gg 1$ limit.

Using Eq.(36) to eliminate the $\bar{x}$ and $\bar{z}$ from Eq.(37) and
Eq.(38) (see Appendix E), we find the coupled equilibrium conditions

$\displaystyle\sum_{j=1}^{D}\frac{1}{\gamma_{j}}\Big{(}\langle
k_{j}\rangle-\langle m_{j}\rangle(1-\delta_{j}^{-1})\Big{)}$
$\displaystyle=\Big{(}N-\langle k\rangle\Big{)}^{2}$ (39)
$\displaystyle\sum_{j=1}^{D}\langle m_{j}\rangle\delta_{j}^{-1}$
$\displaystyle=\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}},$ (40)

where $N\equiv\sum_{j=1}^{D}n_{j}=\sum_{j=1}^{D}r_{j}=D$. Eq.(39) and Eq.(40)
define how the average number of bound and optimally-bound particles for each
species $j$ vary with one another and with the parameters for binding affinity
$\gamma_{j}$ and optimal-binding affinity $\delta_{j}$. For practical
purposes, when trying to solve this system of equations it is necessary to
first solve Eq.(36) and then insert the obtained values of $\bar{z}$ and
$\bar{x}$ into Eq.(37) and Eq.(38) to find $\langle k_{j}\rangle$ and $\langle
m_{j}\rangle$. But Eq.(39) and Eq.(40) do provide an affirming pathway to more
familiar results. If we take $\delta_{j}=\delta$ and $\gamma_{j}=\gamma$ for
all $j$, and note that $\langle k\rangle=\sum_{j=1}^{D}\langle k_{j}\rangle$
(and similarly for $\langle m_{j}\rangle$), we find

$\displaystyle\frac{1}{\gamma}\Big{(}\langle k\rangle-\langle
m\rangle(1-\delta^{-1})\Big{)}=\Big{(}N-\langle
k\rangle\Big{)}^{2},\qquad\langle m\rangle\delta^{-1}=\frac{\langle
k\rangle-\langle m\rangle(1-\delta^{-1})}{N-\langle m\rangle(1-\delta^{-1})}.$
(41)

The results in Eq.(41) are the very same ones we found in [Wil19] for the
gendered dimer system.

By imposing the condition $\langle k\rangle=\langle m\rangle$ on the second
equation in Eq.(41), we can show

$\langle k\rangle=\langle m\rangle=\frac{N-1}{1-\delta^{-1}},$ (42)

suggesting that the condition $\langle k\rangle=\langle m\rangle$ only occurs
when essentially all the particles are bound to their optimal binding sites.
Inserting this value for $\langle k\rangle$ and $\langle m\rangle$ into the
second equation of Eq.(41) yields the thermal condition under which this fully
optimal binding configuration occurs. We find

$N-1=\gamma\delta\frac{(1-N\delta^{-1})^{2}}{1-\delta^{-1}}.$ (43)

In [Wil19], we used Eq.(43) to infer the existence of generally two types of
binding systems with quite different relationships between $\langle k\rangle$
and $\langle m\rangle$. When $\delta\gg\gamma>1$, Eq.(43) became
$\gamma\delta\simeq N$ and we had a ”search-limited” system in which optimal
binding was primarily limited by the ability of particles to find their
optimal binding site in the surrounding volume; when $\gamma\gg\delta>1$,
Eq.(43) became $\delta\simeq N$ and we had a ”combinatorics-limited” system in
which optimal binding was primarily limited by the ability of particles to
avoid the combinatorial sea of suboptimal contacts.

As we increased the temperature in search-limited systems, the value of
$\langle m\rangle$ remained close to the value of $\langle k\rangle$ thus
indicating that such systems could have partial binding to sites but with all
such bindings being optimal. Conversely, in combinatorics-limited systems,
increasing the temperature led to the value of $\langle m\rangle$ being much
lower than the value of $\langle k\rangle$ indicating that when particles were
bound, such binding was likely suboptimal. With some heuristic arguments, we
suggested that biophysical systems are more likely to be of the search-limited
type, but such an inference was limited by the simplicity of our model.

In this work, we want to extend the analysis in this simpler case to one where
there are multiple particle types of various copy number and various binding
and optimal-binding affinities. For this general case, the objective is to
find a condition akin to Eq.(43) that will allow us to distinguish various
binding behaviors in the system and thus tell us if our previous
combinatorics-limited and search-limited framings still apply. Due to its
incorporation of multiple-copy number and type-dependent binding affinities,
this more general case will be more biophysically relevant and could thus
serve as a firmer basis for categorizing biophysical systems as one of the two
types.

But before we consider this most general case, we consider two more specific
cases to build the intuition and techniques for how binding and combinatorics
affect ligand-receptor systems.

## 4 Large $N$ Limits of Special and General Cases

In studying the system modeled by Eq.(31), we will first work through two
special cases that establish the intuition and methods we will later apply to
the most general case. The first special case is that of $\delta_{i}=1$ for
all $i$. This is the case where different ligand-types may have different
binding affinities to the set of receptors, but all receptors are equivalent
from the perspective of a single ligand-type. The second special case is that
of $\gamma_{i}\to\infty$ corresponding to a system where ligands can only
exist as attached to a receptor and where the various microstates consist of
derangements of the ligands amongst the set of receptors. With these two cases
established, we will then consider the general case with no prior assumptions
on the values of $\delta_{i}$ and $\gamma_{i}$.

(a)

(b)

Figure 4: Two limiting cases for the general model. In (a), we have the case
of $\delta_{i}=1$. For this case, all receptors are equivalent and a distinct
ligand species (denoted with a distinct color or border in the figure) of type
$i$ has a binding affinity of $\gamma_{i}$. In (b), we have the case where
$\gamma_{i}\to\infty$. For this case, all ligands are bound to the lattice,
and the system microstates consist of the various ways to permute a list with
repeated elements. The top figure shows the ”correct microstate” (where each
ligand is in its optimal binding site), and the bottom figure shows a
”deranged microstate” (where each ligand is in a suboptimal binding site). The
most general case of ligand-receptor binding is a combination of these two
cases.

###  4.1 Simple Binding Model

One simplification of the most general scenario associated with Eq.(25) is to
have each ligand type have the same binding affinity regardless of to which
receptor it binds. This simplification amounts to taking $\delta_{i}=1$ (or
$\Delta_{i}=0$ by Eq.(21)). Phrased differently, this condition implies that,
for a single type, the optimally and suboptimally bound ligand partition
functions are equal, and thus there is no thermal advantage for a ligand to be
bound to any particular receptor; from the perspective a single type of
ligand, all receptors are thermodynamically identical. However, the ligands of
different types are not thermodynamically identical to each other: Since
$\gamma_{i}$ is not presumed to the be the same for all $i$, each ligand type
has a different binding affinity to an arbitrary receptor. Thus, our system
contains many distinguishable ligands in the presence of many identical
receptors.

The thermal implications of this condition can be found by computing the
partition function. Taking $\delta_{i}\to 1$ in Eq.(24), we can define the new
partition function

$W_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\gamma})\equiv\lim_{\\{\delta_{i}\to
1\\}}{\cal
Z_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta},\boldsymbol{\gamma})}.$
(44)

Using the limit identity $\lim_{\lambda\to
0}\,\lambda^{n}L_{n}^{(\alpha)}\left(\frac{x}{\lambda}\right)=(-1)^{n}{x^{n}}/{n!}$,
the definition Eq.(44) ultimately gives us

$W_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\gamma})=\frac{N_{R}!}{\prod_{j}n_{j}!}\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z^{N_{R}+1}}\,e^{z}\prod_{j=1}^{D}\left(1+\gamma_{j}z\right)^{n_{j}},$
(45)

where we used the definition $N_{R}\equiv\sum_{j=1}^{D}r_{j}$. The factor
${N_{R}!}/{\prod_{j}n_{j}!}$ is an important normalization that ensures that
we have the correct counting of microstates. For example, the $\gamma_{j}\gg
1$ (i.e., mostly bound-ligands) limit should yield (for
$\sum_{j}r_{j}\geq\sum_{j}n_{j}$) a partition function that is proportional to
the number of ways to arrange $n_{i}$ identical objects of type $i$ for
$i=1,\ldots,D$ amongst $r_{1}+\cdots+r_{D}$ sites.

We could have obtained Eq.(45) without use of the limiting case by starting
from the expression

$W_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\gamma})=\sum_{\boldsymbol{k}}I_{\boldsymbol{r},\boldsymbol{k}}\prod_{j=1}^{D}\frac{1}{(n_{j}-k_{j})!}\gamma_{j}^{k_{j}}=\frac{N_{R}!}{\prod_{j}n_{j}!}\sum_{\boldsymbol{k}}\frac{1}{(N_{R}-\sum_{j}k_{j})!}\prod_{j=1}^{D}\binom{n_{j}}{k_{j}}\gamma_{j}^{k_{j}},$
(46)

where the sum for each $k_{j}$ runs from $0$ to $n_{j}$. In the summand of the
first equality of Eq.(46), $I_{\boldsymbol{r},\boldsymbol{k}}$ (defined in
Eq.(11)) is the number of ways to arrange $k_{i}$ bound ligands of type $i$,
for $i=1,\ldots,D$ amongst $r_{1}+\cdots+r_{D}$ receptor sites; the quantity
$1/(n_{j}-k_{j})!$ is the the free-particle partition function for type $j$
ligands (there is no additional factor in the numerator due to Eq.(23)); and
$\gamma_{j}^{k_{j}}$ is the bound-particle partition function for all ligands
of type $j$ that are bound to receptors. By using the contour integral
definition of the inverse factorial, we can show that Eq.(46) is equivalent to
Eq.(45).

Using Eq.(27) and Eq.(46), we find that the average number of bound ligands in
the system can be written as

$\langle
k\rangle=N_{R}\sum_{i=1}^{D}\gamma_{i}\frac{W_{\boldsymbol{n}_{i},\boldsymbol{r}_{i}}(\boldsymbol{\gamma})}{W_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\gamma})},$
(47)

where $\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})$. The
vector $\boldsymbol{r}_{j}$ is defined similarly. Eq.(47) provides us with a
reliable means for computing the order parameter presuming we have a reliable
means for computing the partition function. For our use case, we will use the
large $N$ saddle-point approximation as the basis for this latter computation.

First, defining

${\cal A}_{\boldsymbol{n},\boldsymbol{r}}(z;\boldsymbol{\gamma})\equiv
z-N_{R}\ln
z+\sum_{j=1}^{D}n_{j}\ln\left(1+\gamma_{j}z\right)+\ln\frac{N_{R}!}{\prod_{j}n_{j}!},$
(48)

and then applying the saddle point approximation to
$W_{\boldsymbol{n},\boldsymbol{r}}$ defined in Eq.(45), we obtain the
approximate partition function

$\displaystyle W_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\gamma})$
$\displaystyle=\frac{1}{2\pi i}\oint_{\Gamma}\frac{dz}{z}\exp\big{[}{\cal
A}_{\boldsymbol{n},\boldsymbol{r}}(z;\boldsymbol{\gamma})\big{]}\simeq\frac{1}{\left(2\pi\bar{z}^{2}A^{\prime\prime}_{\boldsymbol{n},\boldsymbol{r}}(\bar{z};\boldsymbol{\gamma})\right)^{1/2}}\exp\left[A_{\boldsymbol{n},\boldsymbol{r}}(\bar{z};\boldsymbol{\gamma})\right],$
(49)

where $\bar{z}$ is defined by the constraint
$0=\partial_{z}A_{\boldsymbol{n},\boldsymbol{r}}(z;\boldsymbol{\gamma})|_{z=\bar{z}}\equiv
A_{\boldsymbol{n}}^{\prime}(z;\boldsymbol{\gamma})|_{z=\bar{z}}$. Eq.(49) is
an approximation, but we henceforth use an equality symbol for notational
sparsity. Using Eq.(48) to compute $\bar{z}$ given its constraint definition,
we find the condition

$0=\bar{z}-N_{R}+\sum_{j=1}^{D}\frac{n_{j}\gamma_{j}\bar{z}}{1+\gamma_{j}\bar{z}}$
(50)

Next, from Eq.(27) and Eq.(44), we can identify the average number of bound
ligands of type $i$, $\langle k_{i}\rangle$, as

$\langle k_{i}\rangle=\gamma_{i}\frac{\partial}{\partial\gamma_{i}}\ln
W_{\boldsymbol{n},\boldsymbol{r}}.$ (51)

From this definition and Eq.(49) and Eq.(48), we thus find

$\langle k_{i}\rangle=\frac{n_{i}\gamma_{i}\bar{z}}{1+\gamma_{i}\bar{z}},$
(52)

where we dropped sub-leading terms. This expression could also have been
derived by starting from Eq.(47) and using Eq.(49) (and dropping the
exponential pre-factor) with the fact that

$A_{\boldsymbol{n}_{i},\boldsymbol{r}_{i}}(\bar{z};\boldsymbol{\gamma})-A_{\boldsymbol{n},\boldsymbol{r}}(\bar{z};\boldsymbol{\gamma})=\ln\bar{z}-\ln(1+\gamma_{i}\bar{z})+\ln\frac{n_{i}}{N_{R}}.$
(53)

Eq.(52) and Eq.(50) provide a means for approximately determining the number
of bound ligands of a specific type. Given the list of binding affinities
$\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{D})$ for our ligands, we can
numerically solve Eq.(50) and insert the result into Eq.(52). However, there
is a different perspective on these results that connects them to the way
binding is typically represented in chemistry. We can use our system of
equations to solve for $\bar{z}$ in terms of $\langle k_{j}\rangle$ in two
ways: We have

$\bar{z}=N_{R}-\langle
k\rangle,\qquad\bar{z}\left(N_{L}-N_{R}+\bar{z}\right)=\sum_{j=1}^{D}\langle
k_{j}\rangle\gamma_{j}^{-1},$ (54)

where we used $N_{L}\equiv\sum_{j=1}^{D}n_{j}$. The first equation is found by
substituting Eq.(52) into Eq.(50). The second equation is found by dividing
Eq.(52) by $\gamma_{i}\bar{z}$, noting that
$n_{j}\gamma_{j}\bar{z}/(1+\gamma_{j}\bar{z})=n_{j}-n_{j}/(1+\gamma_{j}\bar{z})=n_{j}-\langle
k_{j}\rangle/\gamma_{j}\bar{z}$ and inserting this expression into Eq.(50).
Eliminating $\bar{z}$ from Eq.(54), we then have

$\left(N_{R}-\langle k\rangle\right)\left(N_{L}-\langle
k\rangle\right)=\sum_{j=1}^{D}\langle k_{j}\rangle\gamma_{j}^{-1}$ (55)

which is the law of mass action for a system with $N_{R}$ receptors and
$N_{L}$ ligands where the ligand of type $i$ has a binding affinity
$\gamma_{i}$ to any receptor. We could have anticipated this result: When
$\delta_{i}=1$, the combinatorics of various binding configurations becomes
irrelevant since all re-orderings of the same set of bound ligands are
thermally equivalent, and thus the equilibrium properties should be governable
by averages constrained by the law of mass action alone.

The order of our approximation does not allow for $\langle k\rangle=N_{R}$ or
$\langle k\rangle=N_{L}$, which are the cases where all receptors are occupied
or where all ligands are bound, respectively. However, it does allow for
$\langle k\rangle=N_{R}-1$ or $\langle k\rangle=N_{L}-1$ assuming $N_{R}\leq
N_{L}$ and $N_{L}\leq N_{R}$, respectively. These are states where essentially
all receptors are bound with a ligand or all ligands are attached to a
receptor. To move forward, we will assume we have a system where there are
more receptor sites than ligands, i.e., $N_{L}\leq N_{R}$; the alternative
case can be easily analyzed as well. Our almost-completely bound state is then
$\langle k\rangle=N_{L}-1$ and with our definition of $\bar{z}$ in Eq.(54), we
have $\bar{z}=N_{R}-N_{L}+1$. Taking $N_{L}\gg 1$ (as is the case in our large
$N$ approximation), the thermal condition that defines this almost-completely
bound state is therefore

$N_{R}-N_{L}+1=\sum_{j=1}^{D}n_{j}\gamma_{j}^{-1}+O(\gamma_{j}^{-2}).$ (56)

To make this problem easier to analyze, we will constrain ourselves to work
within the case where $n_{j}=r_{j}$ for all $j$, namely matched population of
ligands and receptors. This case is also the one we will explore in the
simulation comparison of this result. For this case, Eq.(56) becomes

$1=\sum_{j=1}^{D}n_{j}\gamma_{j}^{-1}+O(\gamma_{j}^{-2}).$ (57)

Given a temperature-dependence for $\gamma_{j}$, we can numerically solve
Eq.(57) for the temperature at which essentially all of the ligands are bound
to receptor sites. This case of $n_{j}=r_{j}$ for all $j$ is also the case we
will explore in the simulation comparison of the results in this section.

###  4.2 Derangement-Only Model

Another simplification of the most general scenario associated with Eq.(25) is
to consider it as purely a combinatorial one in which the only available
microstates are those consisting of permutations of ligand positions amongst
the receptor sites. This can only occur if the bound-ligand partition function
is infinitely larger that the corresponding unbound-ligand partition function.
Quantitatively, by Eq.(21), this amounts to taking $\gamma_{i}\to\infty$. In
such a case, all ligands are bound to a receptor (if
$\sum_{j}n_{j}<\sum_{j}r_{j}$) or all receptors are occupied (if
$\sum_{j}r_{j}<\sum_{j}n_{j}$), and thermal fluctuations only lead the ligands
to switching receptors.

For simplicity going forward, we will assume $r_{j}\geq n_{j}$ for all $j$.
This corresponds to the situation where there are more receptor sites than
ligands. Taking the partition function Eq.(25) to the $\gamma_{i}\to\infty$
limit and dividing out the thermodynamic pre-factor
$\prod_{i=1}^{D}\gamma_{i}^{n_{i}}$ representing the bound-ligand partition
functions, we can define

$X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})\equiv\lim_{\\{\gamma_{i}\to\infty\\}}\frac{1}{\prod_{i=1}^{D}\gamma_{i}^{n_{i}}}{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta},\boldsymbol{\gamma}),$
(58)

Computing the limit gives us

$X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})=\frac{1}{(N_{R}-N_{L})!}\int^{\infty}_{0}dx\,e^{-x}x^{N_{R}-N_{L}}\prod_{i=1}^{D}(\delta_{i}-1)^{n_{i}}L_{n_{i}}^{(r_{i}-n_{i})}\left(\frac{x}{1-\delta_{i}}\right),$
(59)

where $L_{n}(x)$ is the $n$th Laguerre polynomial. In defining Eq.(58), we
used the fact that the $\\{\gamma_{i}\\}$ are thermodynamically irrelevant for
a system consisting only of bound ligands and thus their factors can be
divided out of the partition function.

We could have derived Eq.(58) without a limiting case by recognizing that the
various microstates of this system are ”partial derangements” of a list.
Specifically, we could have written $X_{\boldsymbol{n},\boldsymbol{r}}$ as a
summation over these derangements:

$X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})=\sum_{\boldsymbol{m}}B_{\boldsymbol{r}-\boldsymbol{m},\boldsymbol{n}-\boldsymbol{m}}\prod_{j=1}^{D}\binom{r_{j}}{m_{j}}\delta_{j}^{m_{j}},$
(60)

where the summation for each $m_{j}$ runs from $0$ to $n_{j}$, and
$B_{\boldsymbol{n},\boldsymbol{k}}$ is defined in Eq.(6). Deriving Eq.(59)
from Eq.(60) requires the identity Eq.(120) derived in Appendix B.2.

Also, Eq.(58) is an ”elements with repeats” generalization of the permutation
glass considered in [Wil18]. If we take $r_{i}=n_{i}=1$ for all $i$ in the
product in Eq.(58), we find

$X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})\Big{|}_{r_{i}=n_{i}=1}=\int^{\infty}_{0}dx\,e^{-x}\prod_{i=1}^{N}\big{(}\delta_{i}-1+x\big{)},$
(61)

which, with $\delta_{i}=e^{\beta\Delta_{i}}$, is identical to the partition
function derived in [Wil18]333In the original paper, we defined our
microstates in terms of energy penalties rather than energy benefits so the
partition function here differs from the original by a multiplicative
constant.. In that work, we derived necessary but not sufficient conditions
for the system to settle into the ”completely correct” microstate, which in
our case corresponds to all ligands being bound in their optimal receptors.
Here we attempt to derive analogous conditions for this more general case.

Using Eq.(27), Eq.(58), Eq.(59), and the identity Eq.(28), we find that the
average number of ligands bound to their optimal receptors is

$\langle
m\rangle=\sum_{j=1}^{D}r_{j}\delta_{j}\frac{X_{\boldsymbol{n}_{j},\boldsymbol{r}_{j}}(\boldsymbol{\delta})}{X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})},$
(62)

where $\boldsymbol{n}_{j}$ is $\boldsymbol{n}$ with 1 subtracted from the
$j$th component (i.e.,
$\boldsymbol{n}_{j}=(n_{1},\ldots,n_{j}-1,\ldots,n_{D})$) and
$\boldsymbol{r}_{j}$ is defined similarly. Eq.(62) tells us that if we have a
consistent means for computing the partition function
$X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})$, we can calculate the
order parameter with little extra work. For this system, the consistent means
we have for computing the partition function is the large $N$ approximation.

To implement this approximation, we first define

$F_{\boldsymbol{n},\boldsymbol{r}}(x;\boldsymbol{\delta})\equiv
x-(N_{R}-N_{L})\ln
x-\sum_{j=1}^{D}\ln\left[\left(\delta_{j}-1\right)^{n_{j}}\,L_{n_{j}}^{(r_{j}-n_{j})}\left(\frac{x}{1-\delta_{j}}\right)\right]+\ln[(N_{R}-N_{L})!]\,.$
(63)

Then applying Laplace’s method to $X_{\boldsymbol{n},\boldsymbol{r}}$ defined
in Eq.(59), we have the approximation

$\displaystyle X_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta})$
$\displaystyle=\int^{\infty}_{0}dx\,\exp\left[-F_{\boldsymbol{n},\boldsymbol{r}}(x;\\{\delta_{i}\\})\right]\simeq\left(\frac{2\pi}{F^{\prime\prime}_{\boldsymbol{n},\boldsymbol{r}}(\bar{x};\boldsymbol{\delta})}\right)^{1/2}\exp\left[-F_{\boldsymbol{n},\boldsymbol{r}}(\bar{x};\boldsymbol{\delta})\right]$
(64)

where $\bar{x}$ is defined by the constraint
$0=\partial_{x}F_{\boldsymbol{n}}(x;\\{\delta_{i}\\})|_{x=\bar{x}}\equiv
F_{\boldsymbol{n}}^{\prime}(x;\\{\delta_{i}\\})|_{x=\bar{x}}$. Eq.(64) is an
approximation, but we henceforth use an equality symbol for notational
sparsity. Applying the constraint
$0=F_{\boldsymbol{n}}^{\prime}(x;\\{\delta_{i}\\})|_{x=\bar{x}}$ to Eq.(63)
and using the recursive Laguerre identity
$u\,\partial_{u}L_{n}^{(\alpha)}(u)=nL_{n}^{(\alpha)}(u)-(n+\alpha)L_{n-1}^{(\alpha)}(u)$,
we find that $\bar{x}$ can be computed from

$\bar{x}=\sum_{j=1}^{D}r_{j}\Biggl{(}1-\frac{\displaystyle
L_{n_{j}-1}^{(\omega_{j})}\left(\bar{\sigma}_{j}\right)}{\displaystyle
L_{n_{j}}^{(\omega_{j})}\left(\bar{\sigma}_{j}\right)}\Biggr{)};\qquad\bar{\sigma}_{j}\equiv\frac{\bar{x}}{1-\delta_{j}},$
(65)

where we defined $\omega_{j}\equiv r_{j}-n_{j}$. Using Eq.(62) with Eq.(64)
(and the fact that we are in the $N\gg 1$ limit), we see that the average
number of optimal bindings is

$\displaystyle\langle m\rangle$
$\displaystyle=\sum_{j=1}^{D}\frac{r_{j}\delta_{j}}{\delta_{j}-1}\frac{\displaystyle
L_{n_{j}-1}^{(\omega_{j})}\left(\bar{\sigma}_{j}\right)}{\displaystyle
L_{n_{j}}^{(\omega_{j})}\left(\bar{\sigma}_{j}\right)}.$ (66)

With Eq.(66), we can determine whether a microstate consisting of all ligands
bound to their optimal receptors can be achieved in this system. Finding the
condition that makes such a microstate possible would require us to look at
the low temperature behavior of the system, but let’s momentarily go in the
opposite direction.

What is the behavior of Eq.(66) when $T$ goes to $\infty$? First, as
$T\to\infty$, $\delta_{i}$ goes to $1$. Given the definition of the Laguerre
polynomial, we can derive the limit $\lim_{\lambda\to
0}\,\lambda^{m}L^{(\alpha)}_{m}\left({x}/{\lambda}\right)=(-1)^{m}{x^{m}}/{m!}.$
Using this limit, we find for $T\to\infty$ that $\langle
m\rangle\to\sum_{j=1}^{D}n_{j}r_{j}/\bar{x}$ and
$\bar{x}\to\sum_{j=1}^{D}r_{j}$. Therefore,

$\lim_{T\to\infty}\langle
m\rangle=\frac{\sum_{j=1}^{D}n_{j}r_{j}}{\sum_{j=1}^{D}r_{j}}$ (67)

There are three things to note about this result: First, the fact that it is
non-zero; second, the vanishing scaling with $r_{j}$ as $r_{j}\to\infty$;
third, the linear scaling with $n_{j}$.

With regard to the first fact, a naive entropic argument might lead us to
think that there would be no optimal bindings at infinite temperature (i.e., a
physical regime where the energy advantage of optimal bindings is irrelevant)
since the macrostate of completely deranged ligands ($\langle m\rangle\simeq
0$) would supposedly have the largest number of microstates and thus the
largest entropy. However, Eq.(67) suggests that the macrostate for completely
deranged bindings is in fact not entropically favored at infinite temperature,
and thus that such a macrostate takes up less configuration space than
seemingly more ordered and constrained macrostates. Indeed, when you have
various types of ligands each of which occurs in large numbers, then it
becomes more constraining to require no ligand to be in its optimal binding
site than it is to have some ligands be optimally bound.

With regard to the second fact, we see that since $r_{j}$ appears in both the
denominator and numerator of Eq.(67), the infinite temperature limit of
$\langle m\rangle$ remains finite as $r_{j}\to\infty$. This means that the
number of optimal bindings that result from random assortment does not change
as the number of available receptor sites increases. Increasing the number of
receptors doesn’t make such optimal binding more likely if the number of
ligands remains constant.

Conversely, we see that Eq.(67) does increase with increasing $n_{j}$ meaning
that increasing the number of ligands in the system does increase the number
of optimal bindings that can occur. This makes sense because with more ligands
in the system there is a greater chance that one of those existing ligands
will make their way to an optimal lattice site.

Pharmacologically, these latter two results imply that if one could only make
one change in order increase the odds of thermally random optimal ligand-
receptor binding, one should increase the concentration of the ligands in the
system rather than increasing the number of receptor sites available to them.
It is important to note, however, that recalling our initial assumption, the
limiting behavior shown in Eq.(67) is valid only for the case where $r_{j}\geq
n_{j}$.

We can obtain an alternative interpretation of the meaning Eq.(67) by writing
it in terms of a covariance. Defining the bar-average of a parameter $O_{j}$
as $\overline{O}=\sum_{j=1}^{D}O_{j}/D$ and the covariance
$\text{cov}(r,n)=\overline{r\,n}-\overline{r}\overline{n}$ we have

$\lim_{T\to\infty}\langle
m\rangle=\bar{n}\left(1+\frac{\text{cov}(r,n)}{\bar{r}\,\bar{n}}\right).$ (68)

The quantity $\bar{n}$ is the average particle-number across all types of
ligands and $\text{cov}(r,n)$ represents how much $r_{j}$ and $n_{j}$ vary
together. Eq.(68) shows that the more that $r_{j}$ and $n_{j}$ vary together,
the larger the lower limit on the number of optimally bound ligands. With more
concurrent variability in the number of ligands and receptors of each type, it
becomes more likely that at least some ligands, just from random assorting,
will be bound to their optimal receptors.

Now, we consider the opposite temperature limit under the frame of a specific
question: At what temperature are all of the ligands bound to their optimal
receptors? For analytical simplicity, we will subsequently work within the
case $n_{j}=r_{j}$ for all $j$, namely matched population of ligands and
receptors. Thus all ligands are bound to their optimal receptors when $\langle
m\rangle=N_{R}=N_{L}$. This case allows us to study the combinatorial
properties of this model more simply without having to consider the various
ways subsets of receptors are occupied.

To answer this question, we first use Eq.(65) and Eq.(66) to obtain the
identity

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

where we defined $N\equiv\sum_{j}^{D}r_{j}=\sum_{j}^{D}n_{j}$ for this case,
and the $\langle m_{j}\rangle$s are the elements of the sum in Eq.(66). When
all ligands are bound to their optimal receptors, we have $\langle
m_{j}\rangle=n_{j}=r_{j}$. Thus, at this desired critical temperature, we have
the condition

$\bar{x}=\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}.$ (70)

To move forward, we will make two assumptions whose consistency we will check
at the end of the calculation. First, we assume that the desired temperature
is sufficiently low that $\delta_{i}\gg 1$ for all $i$. Second, we assume that
$\delta_{i}\gg\bar{x}$. With these assumptions, we find
$L^{(0)}_{n_{j}-1}(\bar{\sigma}_{j})/L^{(0)}_{n_{j}}(\bar{\sigma}_{j})=1-\bar{x}\delta_{j}^{-1}+O(\delta_{j}^{-2})$.
Expanding Eq.(66) to first order in $\delta_{j}^{-1}$ yields

$\langle
m\rangle=N+(1-\bar{x})\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}+O(\delta^{-2}),$ (71)

where $O(\delta^{-2})\equiv\sum_{j=1}^{D}O(\delta_{j}^{-2})$ and
$O(\delta_{j}^{-2})$ represents terms of order $\delta_{j}^{-2}$. At the
temperature at which all ligands are optimally bound, we have $\langle
m\rangle=N$. Thus, Eq.(71) implies $\bar{x}=1+O(\delta_{j}^{-2})$ and by
Eq.(70) we obtain the final condition

$1=\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}+O(\delta^{-2}).$ (72)

Given the temperature dependence for $\delta_{j}$, we can numerically solve
Eq.(72) for the temperature at which all of the ligands are optimally bound to
receptor sites. We will do so in Sec. 5.2 when we simulate this system. To
check consistency with our two initial assumptions (i.e., $\delta_{j}\gg 1$
and $\delta_{j}\gg\bar{x}$), we note that, for the large particle-number
limit, Eq.(72) implies $\delta_{j}>n_{j}\gg 1=\bar{x}+O(\delta_{j}^{-2})$, as
we assumed.

###  4.3 General Case

Having explored various limiting cases, we are now ready for the full case. In
this section, our objective is two fold: First, determine the equations for
both the average number of bound and optimally-bound ligands of each type;
second, use these equations to determine the thermal conditions that define
the system settling into the microstate in which each ligand is bound to its
optimal receptor (i.e., the fully optimally bound). To get to either
objective, we first need to approximate the partition function and compute the
standard observables (Eq.(27)) according to this approximation. We will apply
methods similar to those applied to the limiting cases to analyze this general
case.

Applying the saddle point approximation to Eq.(25), we have

$\displaystyle{\cal
Z}_{\boldsymbol{n},\boldsymbol{r}}(\boldsymbol{\delta},\boldsymbol{\gamma})$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z}\int^{\infty}_{0}dx\,\exp\big{[}{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\big{]}\simeq\frac{1}{({\bar{z}}^{2}\det
H)^{1/2}}\exp\big{[}{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(\bar{z},\bar{x};\boldsymbol{\delta},\boldsymbol{\gamma})\big{]}$
(73)

where

$\displaystyle{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(\bar{z},\bar{x};\boldsymbol{\delta},\boldsymbol{\gamma})$
$\displaystyle\equiv\bar{z}-\bar{x}+(N_{R}-N_{L})\ln(x/z)+\sum_{j=1}^{D}\ln\left[(\gamma_{j}(\delta_{j}-1))^{n_{j}}L_{n_{j}}^{(r_{j}-n_{j})}\left(\bar{\phi}_{j}\right)\right];$
(74) $\displaystyle\quad\bar{\phi}_{j}$
$\displaystyle\equiv\frac{\bar{x}}{1-\delta_{j}}\left(1+\frac{1}{\bar{z}\gamma_{j}}\right).$
(75)

The quantity $H$ is the complex hessian matrix of ${\cal F}_{\boldsymbol{n}}$

$H\equiv\partial_{\alpha,\beta}{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\Big{|}_{x,z=\bar{x},\bar{z}},\qquad$
(76)

where the variables $\alpha$ and $\beta$ can be $x$ or $z$. The critical
points $\bar{z}$ and $\bar{x}$ are defined by the conditions

$0=\partial_{z}{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\Big{|}_{x,z=\bar{x},\bar{z}},\qquad
0=\partial_{x}{\cal
F}_{\boldsymbol{n},\boldsymbol{r}}(z,x;\boldsymbol{\delta},\boldsymbol{\gamma})\Big{|}_{x,z=\bar{x},\bar{z}}.$
(77)

Applying the critical point conditions to Eq.(74), we find

$\bar{z}=N_{R}-N_{L}+\sum_{j=1}^{D}\frac{1}{\bar{z}\gamma_{j}+1}\left(n_{j}-r_{j}\frac{\displaystyle
L_{n_{j}-1}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}\right),\qquad\bar{x}=\sum_{j=1}^{D}n_{j}\left(1-\frac{\displaystyle
L_{n_{j}-1}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}\right),$ (78)

where we defined $\omega_{j}\equiv r_{j}-n_{j}$. The associated values of
$\langle k\rangle$ and $\langle m\rangle$ can then be found by applying
Eq.(27) to the approximated partition function Eq.(73) while ignoring the
subleading pre-factor. Doing so, we obtain

$\langle
k_{j}\rangle=\frac{1}{\bar{z}\gamma_{j}+1}\left(n_{j}\bar{z}\gamma_{j}+r_{j}\frac{\displaystyle
L_{n_{j}-1}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}\right),\qquad\langle
m_{j}\rangle=\frac{r_{j}\delta_{j}}{\delta_{j}-1}\frac{\displaystyle
L_{n_{j}-1}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}{\displaystyle
L_{n_{j}}^{(\omega_{j})}\left(\bar{\phi}_{j}\right)}.$ (79)

With Eq.(79), we can use the solutions for $\bar{x}$ and $\bar{z}$ determined
from Eq.(78) to find the average number of bound and correctly bound ligands
of each type, thus fulfilling the first objective.

For the second objective of determining the thermal conditions for fully
optimal ligand-receptor binding, we first use Eq.(78) and Eq.(79) together to
obtain two equations relating the four quantities:

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

where we took $\langle k\rangle=\sum_{j}\langle k_{j}\rangle$. In Sec. 3.1, we
showed that in the gendered dimer assembly system, the correct assembly
condition yielded the result $\langle k\rangle=\langle
m\rangle=(N-1)/(1-\delta^{-1})=N-1+O(\delta^{-1})$. With this result, we were
then able to find the thermal condition that defined the fully optimized
state. We want to do something similar for the more general case considered in
this section.

We will work within the case $n_{j}=r_{j}$ for all $j$, namely matched
population of ligands and receptors. This case is the most amenable to an
analysis of the competing influences of disorder and binding since we would
not need to consider the multiple ways to select subsets of receptors or
ligands for binding.

We start by defining the state of fully optimal ligand-receptor binding in a
way analogous to the definition in Eq.(42): We assert that the system is in
the state where all ligands are optimally bound to receptors when $\langle
k\rangle$ and $\langle m\rangle$ satisfy

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

where we defined $N\equiv N_{R}=N_{L}$,
$O(\delta^{-1})\equiv\sum_{j=1}^{D}O(\delta_{j}^{-1})$ and
$O(\delta_{j}^{-1})$ represents terms of order $\delta_{j}^{-1}$. To find the
thermal condition that defines fully optimal binding, we need to find the
thermal condition that is consistent with Eq.(81). Applying Eq.(81) to Eq.(80)
and Eq.(79), we find

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

where $O(\delta^{-2})\equiv\sum_{j=1}^{D}O(\delta_{j}^{-2})$ (see Appendix F).
Eq.(82) is the general thermal condition that defines fully optimal binding.
Given the temperature dependences of $\delta_{j}$ and $\gamma_{j}$, we can use
Eq.(82) to determine the temperature at which all ligands of all types settle
into their optimal receptors.

Practically, if we wanted to compute the number of bound and optimally bound
ligands of type $j$ (for the case where $r_{j}=n_{j}$ for all $j$), we would
use the first and second equations in Eq.(79), respectively, assuming our
system is in the $n_{j}\gg 1$ regime. Furthermore, if our system satisfied
$\delta_{j}\gg 1$ at low temperature, then we could use Eq.(82) to compute the
thermal conditions under which the system achieves fully optimal binding. We
perform these calculations for a concrete example in the next section.

## 5 Simulation Comparison

Consider a two-dimensional grid of square lattice sites each of which can be
filled with various colored squares. The colored squares represent the ligands
of the system with a specific color defining a ligand type, and the lattice
sites are the receptors. Each color-type binds optimally to a particular
collection of lattice sites. For pictorial convenience we can arrange the
optimal lattice sites for each particle type such that a figure is created.
This way it is obvious whether our system is in the fully optimally bound
configuration. We depict this system in Fig. 5.

As a clarifying point, the model we developed for ligands binding to receptors
applies equally well to a one-dimensional chain as to an $n$-dimensional grid
as long as both are finite. This is because coordinates on a finite grid can
be mapped one-to-one to a finite list such that a fixed collection of objects
exploring various positions in the multi-dimensional grid is equivalent to
those objects being placed in various orderings in a list.

Figure 5: Depiction of particle binding on a grid. The system is characterized
by different particles (i.e., the colored-squares) each of which occurs in
multiples copies and is associated with a set of ”optimally bound positions”
on the grid. In the figure, the optimal positions for a colored square are the
positions that would lead it to reproduce the microstate on the right. On the
left, some squares are on the grid (i.e., bound) and other squares are off the
grid (i.e., free). In the language of the model, we say that particle species
$i$ has a binding affinity of $\gamma_{i}$ to the grid and an additional
binding affinity factor of $\delta_{i}$ when it is bound to its optimal site.

In what follows, we use this graphical lattice model to present simulation
results for the two limiting cases and the general case discussed in Sec. 4.
For this grid-assembly system, the number of ligands and the corresponding
number of optimal receptors for each type are equal, and thus we study our
system for the case where $n_{i}=r_{i}$.

###  5.1 Simple Binding Model Simulation $(\delta_{i}=1)$

In this section, we affirm the theoretical results in Sec. 4.1 by simulating a
simplification of the grid system in Fig. 5 at various temperatures.

(a) (b)

Figure 6: Grid snapshots and theory vs. simulation for simple ligand-receptor
binding from Sec. 4.1: We defined $\gamma_{j}=(\beta E_{V})^{3/2}e^{\beta
E_{j}}$ where $E_{V}=10^{-3}$ and $E_{j}$ was sampled from a normal
distribution ${\cal N}(\mu_{E},\sigma_{E}^{2})$ with $\mu_{E}=6.0$ and
$\sigma_{E}=2.0$. All energy parameters were taken to be dimensionless. In (a)
we see snapshots of the simulated system at various equilibrium temperatures.
Particles not bound to the grid (i.e., free particles) are not shown. In (b)
we plot a theory vs. simulation comparison for $\langle k\rangle$ as a
function of temperature, and we mark the binding temperature
$k_{B}T_{\text{bind}}$ computed from Eq.(83). The ”Large $N$” values of
$\langle k\rangle$ were computed from Eq.(52) summed over $i$. The
”Simulation” values of $\langle k\rangle$ were computed from the results of a
Metropolis Hastings algorithm where each point is the average of the result of
five simulations. We see that Eq.(83) indeed identifies the temperature at
which the system is nearly completely bound and that Eq.(52) accurately models
the simulation. (See Supplementary Code in Sec. 12 for link to code repository
used to produce this figure.)

The simplification is to assume particles have no ”optimal” position on the
grid (i.e., $\delta_{j}=1$) and thus a single particle has the same binding
affinity to every site on the grid. The system was simulated using the
Metropolis Hastings algorithm in which the microstates transitioned into one
another dependent on free energy differences of the form $\beta
E=\sum_{j}k_{j}\ln\gamma_{j}$, where $k_{j}$ is the number of bound particles
of type $j$ in the microstate and $\gamma_{j}$ is the associated binding
affinity. We allowed for two types of microstate transitions: particle binding
to the grid and particle dissociation from the grid. As is typical for
Metropolis Hastings algorithms, to fully determine the transition
probabilities we also had to incorporate the difference in probabilities of
selecting the particles for the forward- and reverse-transitions between
microstates (See Appendix A for a description of a more general simulation
system and Supplementary Code in Sec. 12 for associated code).

To incorporate an explicit temperature dependence into the system, we set
$\gamma_{j}=(\beta E_{V})^{3/2}e^{\beta E_{j}}$ where $\beta=1/k_{B}T$. The
quantity $E_{V}$ represents the volume-based energy of a free particle in the
system (e.g., $E_{V}\equiv h^{2}/2\pi mV^{2/3}$ for an ideal gas particle of
mass $m$) and thus $(\beta E_{V})^{3/2}$ represents the ratio between the
kinetic partition functions for free and bound particles. The quantity $E_{j}$
represents the binding energy of a particle of type $j$. For simplicity, we
did not assume a $j$ dependence for $E_{V}$.

Taking $k_{B}T_{\text{bind}}=\beta_{\text{bind}}^{-1}$ to be the critical
temperature at which the complete (or, more precisely, ”almost-complete”)
binding state is achieved, Eq.(57) thus became

$1=(\beta
E_{V})^{-3/2}\sum_{j=1}^{D}n_{j}e^{-\beta_{\text{bind}}E_{j}}+O(e^{-2\beta_{\text{bind}}E_{j}}),$
(83)

Solving Eq.(83) for $k_{B}T_{\text{bind}}$ gives us the temperature at which
the thermal advantage of each particle binding to the grid (at any site) is
large enough to overcome the entropic disadvantage of the particle existing
statically in the grid rather than freely in the volume. In a sense, Eq.(83)
defines the thermal condition under which all particles are able to search for
and successfully find the grid in the space they occupy. This ”searching” is
encoded by the product of the $V$ and $n_{j}$ factors in the equation: As $V$
(defined in $E_{V}$) and $n_{j}$ increase, the volume in which a particle must
search increases, and the number of particles doing the searching increases,
respectively. Both increases make it more difficult for the system to settle
into a state in which all particles are bound: Increasing volume increases the
space in which particles must search for the grid; increasing the number of
particles increases the number of units that need to conduct this search
successfully. Thus increasing either of these values makes achieving the
complete binding state more difficult, unless the temperature is lowered
sufficiently so that the binding energy is strong enough to overcome the
entropic disadvantage of having the particles exist freely. It is only below
$k_{B}T_{\text{bind}}$ that the searching entropy succumbs to the energy
advantage and the system settles into its full binding state. We employ this
spatial search metaphor to distinguish this system from one grounded in a
combinatorial search of possible states. We discuss this latter system in the
next section.

To simulate the system, we chose numerical values for all parameters. For
simplicity, we took all energy parameters in the system to be dimensionless.
The values of $E_{j}$ defining $\gamma_{j}$ were sampled from a Gaussian
distribution ${\cal N}(\mu_{E},\sigma_{E}^{2})$ with mean $\mu_{E}=6.0$ and
variance $\sigma_{E}=2.0$. The value of $E_{V}$ was set to $E_{V}=10^{-3}$.
The values of $n_{j}$ were determined directly from Fig. 5: Inspecting the
count of squares for each of the $D=8$ colors and taking each color to be a
particle type, we have $\boldsymbol{n}=\boldsymbol{r}=(9,9,10,5,7,6,3,51)$.

In Fig. 6(a), we show the simulated grid at various equilibrium temperatures.
The particles are colored squares where particle-type is distinguished by
color. Particles not bound to the grid are not shown. The values of $k_{B}T$
are dimensionless because we are taking the energy parameters of the system to
be dimensionless. In the (i) image of Fig. 6(a), we see that all particles are
bound although they are not in their ”correct positions” as defined by the
fully optimally bound state in Fig. 5. This is of course because, with
$\delta_{i}=1$, there is no thermal advantage to being in such entropically
limited positions. As the temperature increases, fewer particles occupy the
grid which confirms the basic intuition that the system should ”melt” at
higher temperatures.

In Fig. 6(b), we plot the theoretical temperature-dependence of $\langle
k\rangle=\sum_{j=1}^{D}k_{j}$ against the simulated temperature-dependence. We
mark the points in the curve that are associated with the grid depictions in
Fig. 6(a). The temperature computed from Eq.(83) is denoted as
$k_{B}T_{\text{bind}}$. We see excellent agreement between the simulation
results and the theoretical results. Moreover, the predicted temperature
computed from Eq.(83) accords with the results of the simulation. Inspecting
(ii) in Fig. 6(a), as we expect, above the critical temperature computed from
Eq.(83), the grid is no longer completely bound.

###  5.2 Derangement-Only Model Simulation $(\gamma_{i}\to\infty)$

In this section, we affirm the theoretical results in Sec. 4.2 by simulating a
simplification of the grid system in Fig. 5 at various temperatures. The
simplification is to consider the system for the case in which all particles
remained on the grid (i.e., $\gamma_{j}\to\infty$) and where state transitions
are confined to particles exchanging positions within one another. The system
was simulated using the Metropolis Hastings algorithm where microstates
transitioned into one another contingent on free energy differences of the
form $\beta E=\sum_{j}m_{j}\ln\delta_{j}$, where $m_{j}$ is the number of
optimally bound particles of type $j$ in the microstate and $\delta_{j}$ is
the additional binding affinity factor for optimal binding. We allowed for
only one type of transition: single-step permutations of particle positions
(See Appendix A for a description of a more general simulation system and
Supplementary Code in Sec. 12 for associated code).

To incorporate temperature into the system, we returned to our original
expression for $\delta_{j}$ in Eq.(27): $\delta_{j}=e^{\beta\Delta_{j}}$. We
recall that $\Delta_{j}$ is the energy-advantage an optimal binding has over
any other binding for a ligand of type $j$. The associated critical
temperature at which all particles were optimally bound was defined as
$k_{B}T_{\text{derang}}$. Taking
$k_{B}T_{\text{derang}}=\beta_{\text{derang}}^{-1}$, Eq.(72) yields

$1=\sum_{j=1}^{D}n_{j}e^{-\beta_{\text{derang}}\Delta_{j}}+O(e^{-2\beta_{\text{derang}}\Delta}).$
(84)

Solving Eq.(84) for $k_{B}T_{\text{derang}}$ gives us the temperature at which
the thermal advantage of each particle settling into its optimal site is large
enough to overcome the entropic disadvantage of choosing that site in the
space of all other combinatorial possibilities. The influence of combinatorics
is encoded by $\boldsymbol{n}=(n_{1},n_{2},\ldots,n_{D})$: As $n_{j}$
increases, the number of possible combinatorial states in the system increases
and thus it becomes more difficult for a ligand to thermally select the
optimal site in a sea of suboptimal ones, unless the temperature is lowered to
diminish how much the combinatorial entropy influences the free energy. It is
only below $k_{B}T_{\text{derang}}$ that combinatorial entropy succumbs to the
energy-advantage of optimal sites, and the system settles into its fully
optimal state.

To simulate the system, we chose numerical values for all parameters. The
values of $\Delta_{j}$ were sampled from a Gaussian distribution ${\cal
N}(\mu_{\Delta},\sigma_{\Delta}^{2})$ with mean $\mu_{\Delta}=4.0$ and
variance $\sigma_{\Delta}=2.0$; energy parameters were taken to be
dimensionless. The values of $n_{j}$ were determined directly from Fig. 5:
Taking each particle to represent a particle type, we have
$\boldsymbol{n}=\boldsymbol{r}=(9,9,10,5,7,6,3,51)$.

(a) (b)

Figure 7: Grid snapshots and theory vs. simulation for derangement model Sec.
4.2: We used $\delta_{j}=e^{\beta\Delta_{j}}$ where $\Delta_{j}$ was sampled
from a normal distribution ${\cal N}(\mu_{\Delta},\sigma_{\Delta}^{2})$ with
$\mu_{\Delta}=4.0$ and $\sigma_{\Delta}=2.0$. All energy parameters were taken
to be dimensionless. In (a) we see snapshots of the simulated system at
various equilibrium temperatures. As temperature increases, the system becomes
more ”deranged,” meaning particles are less likely to assume their optimal
binding sites. In (b) we plot a theory vs. simulation comparison for $\langle
m\rangle$ as a function of temperature, and we mark the binding temperature
$k_{B}T_{\text{derang}}$ computed from Eq.(84). The ”Exact” values of $\langle
m\rangle$ were computed from Eq.(59) and Eq.(62). The ”Large $N$” values of
$\langle m\rangle$ were computed from Eq.(66). The ”Simulation” values of
$\langle m\rangle$ were computed from the results of a Metropolis Hastings
algorithm where each point is the average of the result of five simulations.
We see that Eq.(84) indeed identifies the temperature at which the system is
in its fully optimally bound state and that both Eq.(62) and Eq.(66)
accurately model the simulation. (See Supplementary Code in 12 for link to
code repository used to produce this figure.)

In Fig. 7(a), we show the simulated grid at various equilibrium temperatures.
The particles are colored squares where particle-type is distinguished by
color. The values of $k_{B}T$ are dimensionless because we are taking the
energy parameters of the system to be dimensionless. In the (i) image of Fig.
7(a), we see that all particles are bound in their ”correct positions” as
defined by the fully optimally bound state in Fig. 5. As the temperature
increases, the particles become increasingly ”deranged” from their correct
positions, though we note that even at high temperatures some particles (in
particular the ones with large $n_{j}$) do maintain many of their correct
positions. This latter result is consistent with the discussion following
Eq.(67).

In Fig. 7(b), we plot the theoretical temperature-dependence of $\langle
m\rangle=\sum_{j=1}^{D}\langle m_{j}\rangle$ against the simulated
temperature-dependence. In particular we compare the simulations to the
”Exact” theoretical prediction defined in Eq.(62), and the ”Large $N$”
theoretical prediction defined in Eq.(66). We mark the points in the curve
that are associated with the grid depictions in Fig. 7(a). The temperature
computed from Eq.(84) is denoted as $k_{B}T_{\text{derang}}$. We see excellent
agreement between the simulation results and the theoretical results.
Moreover, the predicted temperature computed from Eq.(84) accords with the
results of the simulation. Inspecting (ii) in Fig. 6(a), as we expect, beyond
the critical temperature the grid begins to show deranged particle states.

Having explored the two limiting cases of the general model, we are now
prepared for the fully general case. We will proceed as we did in these two
example sections: Starting with a theoretical analysis stemming from an
approximation and then finally simulating our results. The objective is to
obtain a condition similar to Eq.(83) and Eq.(84) that implicitly defines the
temperature at which the fully optimally bound state is achieved.

###  5.3 General Case Simulation

In this section, we simulate the system outlined in Sec. 4.3 for the lattice
grid depicted in Fig. 5.

To incorporate temperature into the system we defined
$\delta_{j}=e^{\beta\Delta_{j}}$ and $\gamma_{j}=(\beta E_{V})^{3/2}e^{\beta
E_{j}}$, where $\Delta_{j}$ is the energy advantage for the particle of type
$j$ binding to its optimal lattice site. The quantity $E_{V}$ represents the
volume-based energy of a free particle in the system (e.g., $E_{V}\equiv
h^{2}/2\pi mV^{2/3}$ for an ideal gas particle of mass $m$) and thus $(\beta
E_{V})^{3/2}$ represents the ratio between the kinetic partition functions for
free and bound particles. The quantity $E_{j}$ is the ”base-binding energy” of
the particle of type $j$ to any position on the lattice. For example, if each
particle of type $j$ had a non-zero binding affinity only when bound to its
optimal site, we would have $E_{j}=0$ and $\Delta_{j}>0$ for all $j$. With
these thermal dependences, and taking
$k_{B}T_{\text{crit}}=\beta_{\text{crit}}^{-1}$ to be the critical temperature
at which the system achieves the fully optimally bound state, Eq.(82) becomes

$1=\sum_{j=1}^{D}n_{j}e^{-\beta_{\text{crit}}\Delta_{j}}\left(1+(\beta_{\text{crit}}E_{V})^{-3/2}e^{-\beta_{\text{crit}}E_{j}}\right)+O(e^{-2\beta_{\text{crit}}\Delta_{j}}).$
(85)

Comparing Eq.(85) to Eq.(83) and Eq.(84), it appears that the first is a
combination of the latter two. Moreover, given our ”search” and
”combinatorics” interpretation of Eq.(83) and Eq.(84) respectively, it appears
that Eq.(85) embodies aspects of both limits contingent on various relative
values of the parameters. In the next section, we will explore these
relationships further.

To simulate the system, we chose numerical values for all parameters. The
values of $\Delta_{j}$ were sampled from a Gaussian distribution ${\cal
N}(\mu_{\Delta},\sigma_{\Delta}^{2})$ with mean $\mu_{\Delta}=3.0$ and
variance $\sigma_{\Delta}=1.0$. The values of $E_{j}$ were sampled from a
Gaussian distribution ${\cal N}(\mu_{E},\sigma_{E}^{2})$ with mean
$\mu_{E}=12.0$ and variance $\sigma_{E}=4.0$. The value of $E_{V}$ was set to
$E_{V}=10^{-3}$. The values of $n_{j}$ were determined directly from Fig. 5:
Specifically, taking each color to represent a particle type, we had
$\boldsymbol{n}=\boldsymbol{r}=(9,9,10,5,7,6,3,51)$

(a) (b)

Figure 8: Grid snapshots and theory vs. simulation for general ligand-receptor
binding from Sec. 4.3: We used $\delta_{j}=e^{\beta\Delta_{j}}$ and
$\gamma_{j}=(\beta E_{V})^{3/2}e^{\beta E_{j}}$ where $\Delta_{j}$ was sampled
from a normal distribution ${\cal N}(\mu_{\Delta},\sigma_{\Delta}^{2})$ with
$\mu_{\Delta}=3.0$ and $\sigma_{\Delta}=1.0$, $E_{V}=10^{-3}$, and $E_{j}$ was
sampled from a normal distribution ${\cal N}(\mu_{E},\sigma_{E}^{2})$ with
$\mu_{E}=14.0$ and $\sigma_{E}=2.0$. All energy parameters were taken to be
dimensionless. In (a), we have snapshots of the simulated system at various
temperatures. Particles not bound to the grid are not shown. In (b) we plot a
theory vs. simulation comparison for both $\langle k\rangle$ and $\langle
m\rangle$ as functions of temperature. The critical temperature computed from
Eq.(85) is denoted by $k_{B}T_{\text{crit}}$. The theoretical values of
$\langle k\rangle$ and $\langle m\rangle$ were computed from Eq.(79) and
summed over $j$. The ”Simulation” values of $\langle k\rangle$ and $\langle
m\rangle$ were computed from the results of a Metropolis Hastings algorithm
where each point is the average of the result of five simulations. We see that
Eq.(85) identifies the temperature at which the system is in the fully
optimally bound state and that the theory curves well match the simulation
results. Thus the validity of the Large-$N$ solutions are affirmed in this
case. Finally, we note that the $\langle m\rangle$ curve falls away faster
than the $\langle k\rangle$ curve which is a consequence of the specific
parameter choices we made. In Sec. 6 we explore how different choices could
lead to different behaviors for $\langle m\rangle$ relative to $\langle
k\rangle$. (See Supplementary Code in 12 for link to code repository used to
produce this figure.)

In Fig. 8, we display the results of simulated thermodynamics for our grid
assembly system for this general case. For such equilibrium simulations, we
can choose whatever state transitions we like as long as detailed balance is
satisfied, that is, as long as the ratios between the forward and reverse
transitions are equivalent to the ratios of the Boltzmann factors between the
final and initial states [Kra06]. Therefore when simulating the equilibrium
behavior of the system, we chose state transitions which led to an an
efficient-in-time exploration of the state space even if such transitions were
unphysical. We included three state-transitions: particle binding (i.e, ligand
association to a receptor), particle unbinding (i.e., ligand dissociation from
a receptor), and binding permutation (i.e., bound ligands switching receptor
sites). The binding permutation transition does not occur in real biomolecular
systems, but it was useful for our simulations since it ensured that the
system did not remain trapped in non-equilibrated states at low temperature.
Theoretically with particle binding and unbinding alone, the system should
always find its true equilibrium eventually, but consistently finding such an
equilibrium on the finite time scales of realistic simulations is difficult.
Thus the principal effect of including the binding permutation transition is
to reduce the time needed to simulate these systems.

In Fig. 8(a), we show the simulated grid at various equilibrium temperatures.
In the (i) image of Fig. 8(a), we see that all particles are bound in their
”optimal positions” as defined by the fully optimally bound state in Fig. 5.
As the temperature increases, fewer particles occupy the grid and fewer of the
particles which occupy the grid are in their optimal binding sites which is
consistent with our intuition that the system should lose both binding and
combinatorial order as the temperature increases.

In Fig. 8(b), we plot the theoretical temperature-dependences of $\langle
k\rangle=\sum_{j=1}^{D}\langle k_{j}\rangle$ and $\langle
m\rangle=\sum_{j=1}^{D}\langle m_{j}\rangle$ against their simulated
temperature-dependences. We mark the points in the curve that are associated
with the grid depictions in Fig. 8(a). The temperature computed from Eq.(85)
(i.e., the temperature at which fully optimal binding is achieved) is denoted
as $k_{B}T_{\text{crit}}$. We see excellent agreement between the simulation
results and the theoretical results. Moreover, the computed critical
temperature is in accordance with the results of the simulations. Inspecting
(ii) in Fig. 8(a) (which is at a temperature above the critical temperature)
the system is no longer in its fully optimally bound configuration, as we
should expect.

The $\langle m\rangle$ and $\langle k\rangle$ curves depicted in Fig. 8(b)
represent only one type of relationship between the temperature dependences of
total binding and optimal binding. In this case, we see that as we heat the
system above the critical temperature, the average number of optimally bound
particles falls more quickly than does the average number of total bound
particles. Therefore, slightly above the critical temperature we have a grid-
image such as that depicted in (ii) of Fig. 8(a): particles are mostly bound
but not all of them are in their optimal stites.

This thermal relationship between $\langle k\rangle$ and $\langle m\rangle$
was pre-determined by our parameter choices for $\Delta_{j}$, $E_{j}$, and
$E_{V}$ as was the temperature computed from Eq.(85). In the discussion
following Eq.(43) for the gendered dimer assembly model, we noted how such
parameter choices could lead us to categorize the extremes of that assembly
system as one of two types. In the next section, we will attempt to do
something similar with this general ligand-receptor system.

## 6 Search and Combinatorics Limited Systems

In [Wil19], we found that systems of dimer assembly could often be
characterized as either search-limited or combinatorics-limited contingent on
the relationships between the binding parameters. Given that the system we are
currently studying is a generalization of a version of gendered dimer assembly
(as shown in Sec. 3.1), we can naturally ask if the ligand-receptor binding
system exhibits similar divisions.

We again narrow our analysis to the case of matched populations of ligands and
receptors. Without this assumption, we would have to analyze separately the
cases where $r_{j}<n_{j}$ and $n_{j}<r_{j}$ with little additional benefit in
conceptual understanding. Consider the theoretical condition defining the
microstate where all ligands are optimally bound to a receptor site (rewritten
here from Eq.(82)):

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

To obtain Eq.(86), we took $\delta_{j}\gg 1$. Therefore, all of our system-
limits will necessarily be defined with the base assumption that the optimal
ligand-receptor binding affinity is large. But even with this base assumption,
we can still find two important limits that give us different approximations
for the critical temperature.

Assume first that $\gamma_{j}\gg 1$ for all $j$ at this critical temperature.
Then $1/\gamma_{j}$ is sub-dominant in the second term in the parentheses of
Eq.(86), and we have the approximation

$1=\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}+O(\gamma^{-1})+O(\delta^{-2}).\qquad\text{[Combinatorics
Limiting Condition]}$ (87)

Recalling that $\gamma_{j}$ is the base-binding affinity of a type $j$ ligand
to any receptor, we know that if $\gamma_{j}\gg 1$ for all $j$, then all
ligands have a sufficiently strong binding to any receptor that ligand-
receptor binding will occur even at low temperature. Thus, whether the system
consists entirely of optimal bindings is primarily determined by whether
$\delta_{j}$ is strong enough to bias such optimal bindings over their
combinatorial disadvantage. We call Eq.(87) a ”combinatorics-limiting
condition” (akin to Eq.(72)) for fully optimal binding since the achievement
of fully optimal binding is limited by the influence of the combinatorics on
the system.

On the other hand, assume that $\gamma_{j}\ll 1$ (with
$\gamma_{j}\delta_{j}\gg 1$) for all $j$ at the critical temperature. Then
$1/\gamma_{j}$ is dominant in the second term in the parentheses for Eq.(82)
and we have the approximation

$1=\sum_{j=1}^{D}n_{j}\delta_{j}^{-1}\gamma_{j}^{-1}+O(\gamma^{0})+O(\delta^{-2}).\qquad\text{[Search
Limiting Condition]}$ (88)

In this case, the product $\delta_{j}\gamma_{j}$ represents the net-binding
affinity for an unbound ligand to not merely bind to any receptor but to
specifically bind to its optimal receptor. A small value of $\gamma_{j}$ means
that unbound ligands are not attracted to suboptimal receptors, and thus
without the additional optimal-binding affinity factor $\delta_{j}$, ligands
would generally not bind at all. In particular, since $\gamma_{j}\ll 1$ but
$\delta_{j}\gamma_{j}\gg 1$ the optimal binding affinity is already strong
enough to overcome the combinatorial disadvantage of such bindings because no
other bindings besides the optimal ones are thermally favored. Thus achieving
the fully optimally bound state is not limited by the influence of
combinatorics. Instead, we call Eq.(88) a ”search-limiting condition” (akin to
Eq.(57)) to highlight the fact that for this case, achieving the fully
optimally bound state is primarily limited by ligands’ abilities to search for
their optimal receptor sites in the volume they occupy.

Given temperature dependences for $\delta_{j}$ and $\gamma_{j}$, we can
compute the temperatures corresponding to the approximations Eq.(87) and
Eq.(88). We term these temperatures, respectively, $T_{\text{comb}}$ and
$T_{\text{search}}$. Using the temperature dependences which yielded Eq.(85),
we find Eq.(87) and Eq.(88) become, respectively,

$\displaystyle 1=\sum_{j=1}^{D}n_{j}e^{-\beta_{\text{comb}}\Delta_{j}},\qquad
1=(\beta_{\text{search}}E_{V})^{-3/2}\sum_{j=1}^{D}n_{j}e^{-\beta_{\text{search}}(\Delta_{j}+E_{j})}$
(89)

where $k_{B}T_{\text{search}}=\beta_{\text{search}}^{-1}$ and similarly for
$k_{B}T_{\text{comb}}$, and we dropped sub-leading terms for notational
simplicity. We can use the equations in Eq.(89) to establish upper bounds on
the true critical temperature $T_{\text{crit}}$: From how the temperature
parameter appears in each equation and comparing each to Eq.(85), it is
straightforward to show

$T_{\text{crit}}<T_{\text{comb}},\,T_{\text{search}}.$ (90)

Now why is it important to define these limiting cases at all and then compute
temperatures from them? Because the relative values of the temperatures are
associated with qualitatively different behaviors for the $\langle k\rangle$
and $\langle m\rangle$ curves, and thus computing these temperatures can
immediately provide us with a sense of how the difference $\langle
k\rangle-\langle m\rangle$ varies with temperature.

(a)

(b)

(c)

Figure 9: Search-limited, indeterminate, and combinatorics-limited systems: We
took $\delta_{j}=e^{\beta\Delta_{j}}$ and $\gamma_{j}=(\beta
E_{V})^{3/2}e^{\beta E_{j}}$ where $\Delta_{j}$ was sampled from a normal
distribution ${\cal N}(\mu_{\Delta},\sigma_{\Delta}^{2})$, $E_{V}=10^{-3}$ and
$E_{j}$ was sampled from a normal distribution ${\cal
N}(\mu_{E},\sigma_{E}^{2})$. In (a) with
$(\mu_{\Delta},\sigma_{\Delta})=(4.75,2.0)$ and
$(\mu_{E},\sigma_{E})=(16.0,3.0)$, we have a combinatorics-limited system; In
(b) with $(\mu_{\Delta},\sigma_{\Delta})=(6.75,2.0)$ and
$(\mu_{E},\sigma_{E})=(10.75,3.0)$, we have a system of indeterminate or in-
between type; In (c) with $(\mu_{\Delta},\sigma_{\Delta})=(7.7501,2.0)$ and
$(\mu_{E},\sigma_{E})=(3.0,1.0)$, we have a search-limited system. In the
combinatorics-limited system there is a significant difference between
$\langle k\rangle$ and $\langle m\rangle$ above the critical temperature,
indicating that such systems can have a significant fraction of suboptimally
bound ligands. In the search-limited system there is little difference between
$\langle k\rangle$ and $\langle m\rangle$, indicating that even when ligands
in such systems are partially bound, the ligands are primarily bound
optimally. (See Supplementary Code in 12 for link to code repository used to
produce this figure.)

We show this in Fig. 9. In each figure, we plot simulation vs theory curves
for $\langle k\rangle$ and $\langle m\rangle$ (akin to that displayed in Fig.
8(b)), for various parameter distributions of $\gamma_{i}$ and $\delta_{i}$.
We see that depending on these distributions, $k_{B}T_{\text{comb}}$ and
$k_{B}T_{\text{search}}$ have different relative values and these relative
values can in turn be used to infer properties of the relationship between
$\langle k\rangle$ and $\langle m\rangle$. Specifically, if
$T_{\text{comb}}<T_{\text{search}}$, then above the critical temperature, the
difference $\langle k\rangle-\langle m\rangle$ grows quickly indicating that a
significant fraction of ligands can be bound suboptimally to the set of
receptors (Fig. 9(a)). Conversely if $T_{\text{comb}}>T_{\text{search}}$, then
above the critical temperature, the difference $\langle k\rangle-\langle
m\rangle$ is small, indicating that even when the system consists of only
partially bound ligands, most of these ligands are attached to their optimal
receptor sites (Fig. 9(c)).

We can also use these temperatures to approximately define when a system is
either search-limited or combinatorics-limited. According to the arguments
used to derive Eq.(89) (and as is affirmed by the results in Fig. 9), a system
is combinatorics-limited when

$\text{Combinatorics-limited system:}\quad T_{\text{crit}}\simeq
T_{\text{comb}}$ (91)

and a system is search-limited when

$\text{Search-limited system:}\quad T_{\text{crit}}\simeq T_{\text{search}}.$
(92)

Thus, similar to what was found for dimer system self-assembly [Wil19], we
have found that we can categorize the ligand-receptor system as constrained by
two extremes: A search-limited extreme and a combinatorics-limited extreme.
Given the results of the special cases considered in Sec. 4.1 and Sec. 4.2, we
could also rename the search-limited condition and combinatorics-limited
condition as ”binding-limited” and ”derangement-limited” conditions,
respectively; the binding-limited condition defining whether particles can go
from free-space to attaching to their preferred site on the grid, and the
derangement-limited condition defining whether the particles attach to their
preferred site relative to other sites.

###  6.1 Hints at Non-Equilibrium Behavior

In simulating the general system in Sec. 5.3 and Sec. 6, we included the non-
physical ”binding permutation transition” in order to efficiently explore the
state space. Including such a transition to model equilibrium behavior is
allowed given the constraints of detailed balance (i.e., any transition is
permitted as long as forward and backward transition ratios equate to
Boltzmann factor ratios), but if we want to understand realistic non-
equilibrium behavior—such as how a system of initially free ligands binds over
time to a collection of receptors—we can only use the physical transitions of
binding and unbinding. Limiting our transition choices in this way reveals
another difference between combinatorics and search-limited systems: When only
using physically realistic state transitions, combinatorics-limited systems
were more likely to get trapped in non-equilibrated metastable states than
were search-limited systems. In particular, when we only allowed for particle
binding and unbinding transitions in systems where $\gamma_{j}\gg 1$ and the
temperature satisfied $T<T_{\text{crit}}$ (thus indicating the system should
satisfy $\langle k\rangle=\langle m\rangle\simeq N-1$), the simulated system
did not always find the ”true equilibrium” of fully optimal binding even if
analytical predictions suggested it should.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10: Simulations in time of the $k_{B}T=0.5$ values of $m/N$ for the
systems in Fig. 9 allowing for unphysical (top-row) and physical (bottom-row)
types of transitions. Each simulation began with all particles unbound from
the grid. The black dashed line at $m/N=1$ represents the expected equilibrium
value that all simulations should approach (as inferred from the plots in Fig.
9). For (a), (b), and (c), we used all the same state transitions used to
simulate Fig. 9, namely particle binding, particle unbinding, and particle
permutation. For (d), (e), and (f), we only used the physically relevant state
transitions of particle binding and unbinding. Plots (a) and (d) have the
combinatorics-limited parameters used for Fig. 9(a). Plots (b) and (e) have
the indeterminate system parameters used for Fig. 9(b). Plots (c) and (f) have
the search-limited parameters used for Fig. 9(c). We see that while all system
types approach their expected equilibrium values when using the unphysical
(but efficient-in-time) particle permutation transition, only the search-
limited system reproduces the equilibrium estimate on the finite time scale of
simulations when only particle binding and unbinding are allowed transitions.
This suggests that, for optimal ligand-receptor binding, only search-limited
systems are able to avoid kinetic traps in real non-equilibrium situations.
(See Supplementary Code in 12 for link to code repository used to produce this
figure.)

In Fig. 10, we again simulated the systems whose equilibrium properties are
depicted in Fig. 9, except in this case we only tracked the evolution of the
number of optimally bound particles, denoted $m$, over the course of the
simulation. In each case, the system started from a state consisting entirely
of free particles, and then it evolved according to its transition properties.
The black dashed line at $m/N=1$ represents the approximate predicted
equilibrium value of $\langle m\rangle/N$ for the associated system type at
$k_{B}T=0.5$ (as obtained from the plots in Fig. 9). The colored continuous
lines are various simulations of the system at the given temperature. For Fig.
10(a), Fig. 10(b), and Fig. 10(c), we used the same state-transitions used to
simulate Fig. 9: Particle binding, particle unbinding, and particle
permutation. For Fig. 10(d), Fig. 10(e), and Fig. 10(f), we only used the
physically relevant state transitions of particle binding and unbinding. Plots
Fig. 10(a) and Fig. 10(d) have the combinatorics-limited parameters used for
Fig. 9(a). Plots Fig. 10(b) and Fig. 10(e) have the indeterminate system
parameters used for Fig. 9(b). Plots Fig. 10(c) and Fig. 10(f) have the
search-limited parameters used for Fig. 9(c). We see that while all system
types approach their ”correct” equilibrium value when using the unphysical
(but efficient-in-time) particle-permutation transition, only the search-
limited system reproduces the equilibrium result when only particle binding
and unbinding are allowed transition steps.

The physical explanation for this behavior is simple. If we begin in a state
of free ligands in a system satisfying $\gamma_{j}\gg 1$, then the ligands
have sufficiently strong binding to all receptors to bind to any one of them
and not necessarily to their optimal receptors. Once such ligands are bound,
it is unlikely they will dissociate from these receptors and then bind to
their true optimal receptor because $\gamma_{j}$ is so large. This situation
exists in distinction to that of a search-limited system where $\gamma_{j}\ll
1$ and thus where dissociation from suboptimal receptors is thermodynamically
feasible and ligands only bind strongly (and largely irreversibly) to their
optimal binding sites.

Therefore, combinatorics-limited and search-limited are not merely convenient
equilibrium distinctions between systems. They also typify distinctions in
realistic non-equilibrium behavior. Theoretically, for infinite time, all
systems should reach their true equilibrium regardless of whatever transitions
they manifest. However, achieving such a true equilibrium in finite time is
made difficult when kinetic traps exist in the system. In contrast to
combinatorics-limited systems, search-limited systems have virtually no
kinetic-traps since strong binding only occurs when receptors bind to their
optimal receptor sites which are associated with the true equilibrium. Thus,
in cases where ligands need to specifically bind to certain receptors in
finite time, the systems will necessarily need to be search-limited to ensure
for rapid optimal binding.

This ”prediction” is necessarily a soft one because we have yet to define at a
parameter-level the distinction between the two principal system types. In
[Wil19] we derived necessary but insufficient conditions typifying the
distinction between combinatorics and search-limited systems for dimer
assembly, but, in the current work, the index-dependence of our biophysical
parameters make the analogous such conditions difficult to derive.

## 7 Biophysical Implications

The theoretical investigations of the previous sections were motivated by the
biophysics of ligand-receptor binding. Now we will consider whether our
results can help us better understand this starting point.

For the system of ligand-receptor binding depicted in Fig. 1, the results
Eq.(79) afford us the ability to predict $\langle k_{j}\rangle$ and $\langle
m_{j}\rangle$ for various ligand species given $\gamma_{j}$, $\delta_{j}$,
$n_{j}$, and $r_{j}$. With Eq.(86) and temperature dependences for the binding
affinities, we can determine the temperature at which a system (that has a
matched population of ligands and receptors) settles into the fully optimally
bound configuration. Such a prediction would first require finding values for
$\gamma_{j}$, $\delta_{j}$, and $n_{j}$. Given that $\gamma_{j}$ and
$\delta_{j}$ are ”effective” model parameters representing, respectively, a
ligand of type $j$’s binding affinity to a non-optimal site and the binding
affinity advantage to an optimal site, these values would have to be
approximated from available data on ligand-receptor interactions. A simple
approach to this approximation would amount to taking $\gamma_{j}$ to be the
average of the binding affinities for a particular ligand $j$ to all receptors
(besides the optimal one) in a system and $\delta_{j}$ to be the additional
binding affinity factor to the ligand’s optimal receptor (i.e.,
$\gamma_{j}\delta_{j}$ would be ligand $j$’s absolute binding affinity to this
optimal receptor). However, ligand-receptor affinities, though useful
theoretical quantities for modeling, are notoriously difficult to calculate in
practice [JAVH20] so obtaining accurate predictions of $k_{B}T_{\text{crit}}$
might be similarly challenging. Moreover, the assumption of matched ligand and
receptor populations, though useful for obtaining intuitive analytical
results, is clearly not a general one.

Still, we can use the theoretical properties of Eq.(86) to make qualitative
statements about the properties of optimally bound systems. First, we
recognize that the terms within the sum must each be less than unity.
Therefore as one of the factors of a single term increases, the other factors
must decrease in order to keep the total product less than unity. For a system
with a large value of $n_{j}$ for ligand $j$, this means we need a
correspondingly large $\delta_{j}$ in order for the thermal constraint
condition to be soluble. Conceptually, this implies that the more copies we
have of a ligand-type, then the more strongly that ligand must bind to its
optimal receptor in order for the entire system’s fully optimally bound
configuration to be achievable.

This result is also true for the total number of ligands in the system: The
larger the number of total ligands in the system, the greater the average
optimal binding affinities of the ligands must be in order for the system to
settle into its fully optimally bound configuration at constant temperature.
We can see this by applying Jensen’s inequality to the first equation in
Eq.(89) and using Eq.(90). Doing so we obtain

$k_{B}T_{\text{crit}}<\frac{1}{\ln N}\sum_{j=1}^{D}\frac{n_{j}}{N}\Delta_{j},$
(93)

where we defined $N\equiv\sum_{j=1}^{D}n_{j}$. Eq.(93) sets an upper limit on
the the temperature at which all ligands settle into their optimal receptors.
Since $n_{j}/N$ is the fraction of elements of type $j$ and $\Delta_{j}$ is
the associated binding energy benefit for being in an optimal site, Eq.(93)
shows that the critical temperature is bounded above by the weighted average
of energy benefits, $\overline{\Delta}\equiv\sum_{j=1}^{D}n_{j}\Delta_{j}/N$.
Moreover, if we want this bound to remain the same as we increase $N$, the
average optimal binding energy must increase in tandem. Thus the more ligands
we have in the system the greater we expect the average binding affinity to
be, presuming the energetically optimal binding configuration is a desired
state in the system.

One could likely have guessed a result of the form in Eq.(93). Namely, it
makes sense that the critical temperature should be of the same order as an
average binding energy in the system. However, what is perhaps surprising is
that the limiting temperature scales as $1/\ln N$: As the total number of
ligands in the system increases, the temperature at which the low energy
system is accessible decreases in tandem but does so logarithmically. Though
seemingly novel, this scaling appears to be archetypal for combinatorial
statistical physics systems where state spaces are combinations and
permutations of a set [Wil18, Wil19].

We could also invert this inequality and use it to establish an upper-limit on
the number of ligands in the system. Assuming we know $k_{B}T_{\text{crit}}$,
we can impose an upper limit on $N$ as

$N<e^{\beta_{\text{crit}}\bar{\Delta}}.$ (94)

As $\bar{\Delta}$ increases, so too does the limit on $N$. Thus systems with
larger average optimal binding energies can also admit more ligands and still
achieve the fully optimally bound configuration at nearly the same
temperatures. The main biophysical implication is that cells with more
proteins should also have larger average binding affinities for those
proteins. For example, human cells and prokaryotic cells, which differ in
numbers of proteins by a few order of magnitude [Mil13], should also differ in
the average optimal binding energies for such proteins. In particular if the
compared collection of cells are typically found in similar thermal
environments (i.e., exist at the same $k_{B}T$) and the particle number was
stringently bound by Eq.(94), we would expect the relative average values of
the binding energies for human and prokaryotic ligands should be the same as
the order of magnitude that differentiates their relative number of proteins.

However, the limiting effect of Eq.(94) is only relevant if the associated
system is specifically combinatorics-limited, i.e., if Eq.(87) well
approximates the general thermal condition Eq.(86). The results in Sec. 6.1
suggest that real biophysical systems (where optimal ligand-receptor binding
is functionally important) are likely search-limited rather than
combinatorics-limited since it is only in the former that fully optimal
binding can be achieved through physical transitions on finite time scales. We
recall that search-limited systems are those for which the energy advantage
for optimal binding is sufficiently high that the only limiting factor to
ligands finding their optimal site is their ability to ”search” for these
sites in the constituent volume. More formally, search-limited systems are
roughly defined as those for which $T_{\text{search}}$ (computed from Eq.(88))
provides a good approximation for $T_{\text{crit}}$ (computed from Eq.(86))

Given that real ligand-receptor systems are likely search-limited, they would
also likely easily satisfy the combinatorics condition Eq.(87) at their
typical temperature. Therefore, the derived inequality Eq.(94) would not be a
strong limit on the number of particles in systems where optimal binding is
functionally required. Instead, a stronger limit would be found by using the
search-limiting condition Eq.(88) to find a bound on particle number. Finding
such a bound first requires us to determine a form for $\gamma_{j}$ that is a
convex function of temperature, and such a $\gamma_{j}$ depends on our exact
model of how free ligands and bound ligands exist in the system of interest.
But as a toy-case, we can also use our grid-assembly expression for
$\gamma_{j}$ to get a sense of the general form of the limit444We note that