Derangement model of ligand-receptor binding

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.


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 1.0 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-tomany 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 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.
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 ligandreceptor 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.

2.0
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.

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]: 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, . . ., and n D elements of type D (where D ≤ N ) is where n ≡ (n 1 , n 2 , . . . , n D ) and L n (x) is the Laguerre polynomial defined as L n (x) = n j=0 n k 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, . . ., 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 ≤ r 1 elements of type 1, k 2 ≤ r 2 elements of type r 2 , . . ., and k D ≤ 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 r,k where r = (r 1 , r 2 , . . . , r D ) and k = (k 1 , k 2 , . . . , 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 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 We can define analogous factors for B r,k and use them to write a summation expression for the 2.0 quantity. The result is To understand Eq.(4) we consider how each factor in the summand contributes to the final expression. The factor ri ji , for i = 1, . . . , 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) j1+···+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 + · · · + r D − j D possible positions which we must fill with k 1 − j 1 + · · · + 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 ! is the number of ways to permute the k 1 − j 1 + · · · + 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: Using Eq.(5) and the definition of the Gamma function, we find that Eq.(4) can be written as where we defined α i ≡ n i − k i .
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 ≡ 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 ≤ D elements to arrange amongst the lattice sites such that none is in its associated site is Thus, we should be able to show that Eq.(6) reduces to Eq.(7) under the right conditions. In particu- 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 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 r−m,k−m . lar if we take r = (1, 1, . . . , 1) ≡ r 0 (i.e., we have D unique elements, each of a single copy-number), then the vector k in Eq.(6) can only have elements of 1 or 0, and thus k defines a particular subset of the total set of elements. B r0,k then represents the number of ways to completely derange a particular collection of unique elements where the collection is defined by the vector 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 r0,k over all possible values of k such that j k j = K. Thus the consistency check we must make is 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, . . ., k D elements of type D to be placed amongst the available r 1 + · · · + 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 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 r,k , we note that r1 m1 · · · r D m D B r−m,k−m is the number of ways to choose m i out of r i positions (for i = 1, . . . , D) to contain their original 3.0 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 we should obtain the number of ways to arrange (and not necessarily derange) k i ≤ r i elements of type i for i = 1, . . . , D across a total of r 1 + . . . + 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 + · · · + r D sites: There are k 1 + · · · + k D filled sites and r 1 − k 1 + · · · + r D − k D empty sites. Consequently there are (r 1 + · · · + 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 + · · · + r D − k D )!, the number of ways to reorder these empty sites. Thus we should find 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.

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 predefined 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, . . . , 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 ≡ r 1 + · · · + r D total receptors and each of the N L ≡ n 1 + · · · + 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 C with elements C i,j defined as C i,j = # of bindings between ligand of type i and receptor of type j.
If we can specify all elements of the D × D matrix 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 C i,j : Each element must be an integer, Next, we ask how we can incorporate binding parameters to represent the way energy affects the probability of a microstate. We will take Q B i,j 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 ) Ci,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 F i to be the single-particle partition function for an unbound ligand of type i, and given that n i − D j=1 C i,j is the total number of unbound ligands of type i, the multi-particle partition function for the ligand of 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 where n ≡ (n 1 , . . . , n D ) and r = (r 1 , . . . , 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 C i,j according to the constraints of the physical system (i.e., C i,j is an integer, D j=1 C i,j ≤ n i and D i=1 C i,j ≤ 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 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 B i,j . 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 3.0 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 ∆ i (with ∆ i ≥ 0). Mathematically, we can encode this assumption into the model by making the transformation where δ i,j is the Kronecker delta, and β = 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 β∆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 B i,j , we have 2N parameters from (Q B i ) and ∆ i together, a reduction that makes our modeling more tractable. Regarding the physical motivation of this assumption, Eq.
Ci,j δi,j From, here we make a change of variables motivated by the quantities that appear in the expression.
We define k i = D j=1 C i,j as the total number of bound (optimally or not) ligands of type i and m i = D j=1 C i,j δ j,i as the total number of optimally bound ligands of type i. Eq.(15) then becomes where we defined the summations as and the combinatorial factor as resulting from the change of summation variables. Qualitatively Ω(k, r, m) is the number of ways to fill m i out of r i receptor sites (for i = 1, . . . , D) with their optimal binding partner ligands while having the remaining k i − m i bound ligands (again for i = 1, . . . , 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., {C i,j }). However, it seems that in order to compute Eq.(16), we would need to compute Ω(k, r, 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 Ω, we can assert that where B n,k is defined in Eq.
We recall that in the summands of Eq.(20), k i ≤ n i represents the total number of ligands of type i that are bound to receptors, and m i ≤ 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 thus giving us 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 3.0 function ratios) should affect the physics of a system. Thus, without loss of generality, we can impose Q F i = 1 for all i under the assumption that the thermal dependence of each Q F i can be absorbed into a redefinition of γ i and δ i with no change in the physical implications of Eq.(22).
With this imposition we have c n = 1 [free-particle partition function normalization] 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 n,k in Eq.(6) and a few Laguerre polynomial identities. After some work (see Appendix C), we obtain We can write this result in a more mathematically useful form by expressing the integrand as the exponential of a potential function: where Γ is a closed contour about the origin in the complex plane and with L   n−1 (u) allows us to prove 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-toone 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].
We can then show Thus from Eq.(27), we have where n j is n with 1 subtracted from the jth component: n j = (n 1 , . . . , n j − 1, . . . , n D ). The vector r j is defined similarly. Although the term r i − 1 does not appear in Eq.(29), we had to introduce r j into the expression Eq.(30) to ensure that r j −n j in Eq.(24) remained unchanged when we replaced n with n j . Eq.(30) could also have been derived from Eq.(22) by differentiating with respect to δ j and making the dummy variable replacement k j = k j − 1 and m j = m j − 1. There is no simplified expression for k 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.

Gendered Dimer-System Assembly
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 δ i and γ 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 where 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): Applying the large N integral approximation (specifically D 1 in this case) to Eq.(31) yields where H is the hessian matrix with second-order derivative components H α,β = ∂ α,β F D z=z,x=x (α, β ∈ {x, z}), andz andx are defined by the conditions In order to compute equilibrium conditions for k j and m j from Eq.(34) we first need the conditions forx andz 1 . Using Eq.(32) and Eq.(35) to find these conditions, we have, from ∂ x F D = 0 and Next, computing k j and m j , we have where, in applying Eq.(33) to Eq.(34), we neglected the exponential pre-factor in the latter since it is subleading in the D 1 limit.
Using Eq.(36) to eliminate thex andz from Eq.(37) and Eq.(38) (see Appendix E), we find the 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 γ j and optimal-binding affinity δ 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 ofz and x into Eq.(37) and Eq.(38) to find k j and m j . But Eq.(39) and Eq.(40) do provide an affirming pathway to more familiar results. If we take δ j = δ and γ j = γ for all j, and note that k = D j=1 k j (and similarly for m j ), we find The results in Eq.(41) are the very same ones we found in [Wil19] for the gendered dimer system.
By imposing the condition k = m on the second equation in Eq.(41), we can show suggesting that the condition k = m only occurs when essentially all the particles are bound to their optimal binding sites. Inserting this value for k and m into the second equation of Eq. (41) yields the thermal condition under which this fully optimal binding configuration occurs. We find In [Wil19], we used Eq.(43) to infer the existence of generally two types of binding systems with quite different relationships between k and m . When δ γ > 1, Eq.(43) became γδ 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 γ δ > 1, Eq.(43) be-

4.1
came δ 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 m remained close to the value of k 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 m being much lower than the value of k 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 typedependent 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.

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 δ 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 γ i → ∞ 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 δ i and γ i .

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 δ i = 1 (or ∆ 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 γ 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 δ i → 1 in Eq.(24), we can define the new partition function Using the limit identity lim λ→0 λ n L (α) n x λ = (−1) n x n /n!, the definition Eq.(44) ultimately gives us where we used the definition N R ≡ D j=1 r j . The factor N R !/ j n j ! is an important normalization that ensures that we have the correct counting of microstates. For example, the γ j 1 (i.e., mostly bound-ligands) limit should yield (for j r j ≥ 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, . . . , D amongst r 1 + · · · + r D sites.
We could have obtained Eq.(45) without use of the limiting case by starting from the expression where the sum for each k j runs from 0 to n j . In the summand of the first equality of Eq.(46), I r,k (defined in Eq.(11)) is the number of ways to arrange k i bound ligands of type i, for i = 1, . . . , D amongst r 1 +· · ·+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 γ kj j is the

Simple Binding Model
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 where n j is n with 1 subtracted from the jth component: n j = (n 1 , . . . , n j − 1, . . . , n D ). The vector 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 and then applying the saddle point approximation to W n,r defined in Eq. (45), we obtain the approximate partition function wherez is defined by the constraint 0 = ∂ z A n,r (z; γ)| z=z ≡ A n (z; γ)| z=z . Eq.(49) is an approximation, but we henceforth use an equality symbol for notational sparsity. Using Eq.(48) to computez given its constraint definition, we find the condition Next, from Eq.(27) and Eq.(44), we can identify the average number of bound ligands of type i, k i , as From this definition and Eq.(49) and Eq.(48), we thus find 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 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 γ = (γ 1 , . . . , γ 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 forz in terms of k j in two ways: We havē where we used N L ≡ D j=1 n j . The first equation is found by substituting Eq.(52) into Eq.(50). The second equation is found by dividing Eq.(52) by γ iz , noting that n j γ jz /(1 + γ jz ) = n j − n j /(1 + γ jz ) = n j − k j /γ jz and inserting this expression into Eq.(50). Eliminatingz from Eq.(54), we then have 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 γ i to any receptor. We could have anticipated this result: When δ i = 1, Taking N L 1 (as is the case in our large N approximation), the thermal condition that defines this almost-completely bound state is therefore 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 Given a temperature-dependence for γ 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.

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 γ i → ∞. In such a case, all ligands are bound to a receptor (if j n j < j r j ) or all receptors are occupied (if j r j < j n j ), and thermal fluctuations only lead the ligands to switching receptors.
For simplicity going forward, we will assume r j ≥ 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 γ i → ∞ limit and dividing out the thermodynamic pre-factor D i=1 γ ni i representing the bound-ligand partition functions, we can define Computing the limit gives us where L n (x) is the nth Laguerre polynomial. In defining Eq.(58), we used the fact that the {γ 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 n,r as a summation over these derangements: where the summation for each m j runs from 0 to n j , and B n,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 which, with δ i = e β∆i , is identical to the partition function derived in [Wil18] 2 . 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 where n j is n with 1 subtracted from the jth component (i.e., n j = (n 1 , . . . , n j − 1, . . . , n D )) and r j is defined similarly. Eq.(62) tells us that if we have a consistent means for computing the partition function X n,r (δ), 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 Then applying Laplace's method to X n,r defined in Eq.(59), we have the approximation wherex is defined by the constraint 64) is an approximation, but we henceforth use an equality symbol for notational sparsity. Applying the constraint 0 = F n (x; {δ i })| x=x to Eq.(63) and using the recursive Laguerre identity u ∂ u L (α) n−1 (u), we find thatx can be computed from where we defined ω j ≡ r j − n j . Using Eq.(62) with Eq.(64) (and the fact that we are in the N 1 limit), we see that the average number of optimal bindings is 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 ∞? First, as T → ∞, δ i goes to 1. Given the definition of the Laguerre polynomial, we can derive the limit lim λ→0 λ m L (α) Using this limit, we find for 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 → ∞; 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 ( m 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 m remains finite as r j → ∞. 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 ≥ n j .
We can obtain an alternative interpretation of the meaning Eq.
The quantityn is the average particle-number across all types of ligands and 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 m = 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 identitȳ where we defined N ≡ D j r j = D j n j for this case, and the m j s are the elements of the sum in Eq.(66). When all ligands are bound to their optimal receptors, we have m j = n j = r j . Thus, at this desired critical temperature, we have the condition 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 δ i 1 for all i. Second, we assume that δ i x. With these assumptions, we find L where O(δ −2 ) ≡ Given the temperature dependence for δ 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., δ j 1 and δ j x), we note that, for the large particle-number limit, Eq.(72) implies δ j > n j 1 =x+O(δ −2 j ), as we assumed.

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 4.3 General Case 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 where The quantity H is the complex hessian matrix of F n where the variables α and β can be x or z. The critical pointsz andx are defined by the conditions 0 = ∂ z F n,r (z, x; δ, γ) x,z=x,z , 0 = ∂ x F n,r (z, x; δ, γ) x,z=x,z .
Applying the critical point conditions to Eq.(74), we find where we defined ω j ≡ r j − n j . The associated values of k and m 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 With Eq.(79), we can use the solutions forx andz 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 quanti- where we took k = j k j . In Sec. 3.1, we showed that in the gendered dimer assembly system, the correct assembly condition yielded the result k = m = (N − 1)/(1 − δ −1 ) = N − 1 + O(δ −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 k and m satisfy 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 where (82) is the general thermal condition that defines fully optimal binding. Given the temperature dependences of δ j and γ 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 1 regime. Furthermore, if our system satisfied δ j 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.

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 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 γ i to the grid and an additional binding affinity factor of δ i when it is bound to its optimal site. objects exploring various positions in the multi-dimensional grid is equivalent to those objects being placed in various orderings in a list.
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 .

Simple Binding Model Simulation (δ 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. The simplification is to assume particles have no "optimal" position on the grid (i.e., δ 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 βE = j k j ln γ j , where k j is the number of bound particles of type j in the microstate and γ 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 γ j = (βE V ) 3/2 e βEj where β = 1/k B T . The quantity E V represents the volume-based energy of a free particle in the We defined γ j = (βE V ) 3/2 e βEj where E V = 10 −3 and E j was sampled from a normal distribution N (µ E , σ 2 E ) with µ E = 6.0 and σ 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 k as a function of temperature, and we mark the binding temperature k B T bind computed from Eq.(83). The "Large N " values of k were computed from Eq.(52) summed over i. The "Simulation" values of k 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.) system (e.g., E V ≡ h 2 /2πmV 2/3 for an ideal gas particle of mass m) and thus (β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 bind = β −1 bind to be the critical temperature at which the complete (or, more precisely, "almost-complete") binding state is achieved, Eq.(57) thus became Solving Eq.(83) for k B T 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

Derangement-Only Model Simulation (γ i → ∞)
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 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 γ j were sampled from a Gaussian distribution N (µ E , σ 2 E ) with mean µ E = 6.0 and variance σ 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 n = r = (9, 9, 10, 5, 7, 6, 3, 51).
In Fig. 6a, 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. 6a, 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 δ 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. 6b, we plot the theoretical temperature-dependence of k = D j=1 k j against the simulated temperature-dependence. We mark the points in the curve that are associated with the grid depictions in Fig. 6a. The temperature computed from Eq.(83) is denoted as k B T 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. 6a, as we expect, above the critical temperature computed from Eq.(83), the grid is no longer completely bound.

Derangement-Only Model Simulation (γ i → ∞)
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., γ j → ∞) 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 βE = j m j ln δ j , where m j is the number of optimally bound particles of type j in the microstate and δ 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 δ j in Eq.(27): δ j = e β∆j . We recall that ∆ 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 derang . Taking k B T derang = β −1 derang , Eq.(72) yields Solving Eq.(84) for k B T 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 n = (n 1 , n 2 , . . . , 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 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 ∆ j were sampled from a Gaussian distribution N (µ ∆ , σ 2 ∆ ) with mean µ ∆ = 4.0 and variance σ ∆ = 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 n = r = (9, 9, 10, 5, 7, 6, 3, 51).
In Fig. 7a, 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. 7a, 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. 7b, we plot the theoretical temperature-dependence of m = D j=1 m j 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. 7a. The temperature computed from Eq.(84) is denoted as k B T 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. 6a, 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. We used δ j = e β∆j where ∆ j was sampled from a normal distribution N (µ ∆ , σ 2 ∆ ) with µ ∆ = 4.0 and σ ∆ = 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 m as a function of temperature, and we mark the binding temperature k B T derang computed from Eq.

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 δ j = e β∆j and γ j = (βE V ) 3/2 e βEj , where ∆ 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 ≡ h 2 /2πmV 2/3 for an ideal gas particle of mass m) and thus (β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 ∆ j > 0 for all j. With these thermal dependences, and taking k B T crit = β −1 crit to be the critical temperature at which the system achieves the fully optimally bound state, Eq. 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 m curve falls away faster than the k 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 m relative to k . (See Supplementary Code in 12 for link to code repository used to produce this figure.) spectively, 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 ∆ j were sampled from a Gaussian distribution N (µ ∆ , σ 2 ∆ ) with mean µ ∆ = 3.0 and variance σ ∆ = 1.0. The values of E j were sampled from a Gaussian distribution N (µ E , σ 2 E ) with mean µ E = 12.0 and variance σ 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 n = r = (9, 9, 10, 5, 7, 6, 3, 51) 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 6.0 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. 8a, we show the simulated grid at various equilibrium temperatures. In the (i) image of Fig. 8a, 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. 8b, we plot the theoretical temperature-dependences of k = D j=1 k j and m = D j=1 m j against their simulated temperature-dependences. We mark the points in the curve that are associated with the grid depictions in Fig. 8a. The temperature computed from Eq.(85) (i.e., the temperature at which fully optimal binding is achieved) is denoted as k B T 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. 8a (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 m and k curves depicted in Fig. 8b 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. 8a: particles are mostly bound but not all of them are in their optimal stites. This thermal relationship between k and m was pre-determined by our parameter choices for ∆ 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.

Search and Combinatorics Limited Systems
In [Wil19], we found that systems of dimer assembly could often be characterized as either searchlimited 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)): To obtain Eq.(86), we took δ j 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 γ j 1 for all j at this critical temperature. Then 1/γ j is sub-dominant in the second term in the parentheses of Eq.(86), and we have the approximation Recalling that γ j is the base-binding affinity of a type j ligand to any receptor, we know that if γ j 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 δ 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 γ j 1 (with γ j δ j 1) for all j at the critical temperature. Then 1/γ j is dominant in the second term in the parentheses for Eq.(82) and we have the approximation In this case, the product δ j γ 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 γ j means that unbound ligands are not attracted to suboptimal receptors, and thus without the additional optimalbinding affinity factor δ j , ligands would generally not bind at all. In particular, since γ j 1 but δ j γ j 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. Figure 9: Search-limited, indeterminate, and combinatorics-limited systems: We took δ j = e β∆j and γ j = (βE V ) 3/2 e βEj where ∆ j was sampled from a normal distribution N (µ ∆ , σ 2 ∆ ), E V = 10 −3 and E j was sampled from a normal distribution N (µ E , σ 2 E ). In (a) with (µ ∆ , σ ∆ ) = (4.75, 2.0) and (µ E , σ E ) = (16.0, 3.0), we have a combinatorics-limited system; In (b) with (µ ∆ , σ ∆ ) = (6.75, 2.0) and (µ E , σ E ) = (10.75, 3.0), we have a system of indeterminate or in-between type; In (c) with (µ ∆ , σ ∆ ) = (7.7501, 2.0) and (µ E , σ E ) = (3.0, 1.0), we have a search-limited system. In the combinatorics-limited system there is a significant difference between k and m 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 k and m , indicating that even when ligands in such systems are partially bound, the ligands are primarily bound optimally.
where k B T search = β −1 search and similarly for k B T 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 crit : From how the temperature parameter appears in each equation and comparing each to Eq.(85), it is straightforward to show 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 k and m curves, and thus computing these temperatures can immediately provide us with a sense of how the difference k − m varies with temperature.
We show this in Fig. 9. In each figure, we plot simulation vs theory curves for k and m (akin to that displayed in Fig. 8b), for various parameter distributions of γ i and δ i . We see that depending on these distributions, k B T comb and k B T search have different relative values and these relative values can in turn be used to infer properties of the relationship between k and m . Specifically, if T comb < T search , then above the critical temperature, the difference k − m grows quickly indicating that a significant fraction of ligands can be bound suboptimally to the set of receptors (Fig. 9a). Conversely if T comb > T search , then above the critical temperature, the difference k − m 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. 9c).
We can also use these temperatures to approximately define when a system is either searchlimited 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 Combinatorics-limited system: T crit T comb (91) and a system is search-limited when Search-limited system: T crit T search .
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 bindinglimited 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.

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 γ j 1 and the temperature satisfied T < T crit (thus indicating the system should satisfy k = m N − 1), the simulated system did not always find the "true equilibrium" of fully optimal binding even if analytical predictions suggested it should.
In Fig. 10, we again simulated the systems whose equilibrium properties are depicted in Fig. 9,  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. 9a. Plots (b) and (e) have the indeterminate system parameters used for Fig. 9b. Plots (c) and (f) have the search-limited parameters used for Fig.  9c. 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.) 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 m /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. 10a, Fig. 10b, and Fig. 10c, we used the same state-transitions used to simulate Fig. 9: Particle binding, particle unbinding, and particle permutation. For Fig. 10d, Fig. 10e, and Fig. 10f, we only used the physically relevant state transitions of particle binding and unbinding. Plots Fig. 10a and Fig. 10d have the combinatoricslimited parameters used for Fig. 9a. Plots Fig. 10b and Fig. 10e have the indeterminate system parameters used for Fig. 9b. Plots Fig. 10c and Fig. 10f have the search-limited parameters used for Fig. 9c. 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.

Hints at Non-Equilibrium Behavior
The physical explanation for this behavior is simple. If we begin in a state of free ligands in a system satisfying γ j 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 γ j is so large. This situation exists in distinction to that of a search-limited system where γ j 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.

Biophysical Implications
The theoretical investigations of the previous sections were motivated by the biophysics of ligandreceptor 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 k j and m j for various ligand species given γ j , δ 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 γ j , δ j , and n j . Given that γ j and δ 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 γ 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 δ j to be the additional binding affinity factor to the ligand's optimal receptor (i.e., γ j δ j would be ligand j's absolute bind-7.0 ing 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 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 δ 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 where we defined N ≡ D j=1 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 ∆ 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, ∆ ≡ D j=1 n j ∆ 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 crit , we can impose an upper limit on N as As∆ 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. 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 search (computed from Eq.(88)) provides a good approximation for T 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 searchlimiting condition Eq.(88) to find a bound on particle number. Finding such a bound first requires us to determine a form for γ j that is a convex function of temperature, and such a γ 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 γ j to get a sense of the general form of the limit 3 . Using the grid-assembly expression for γ j amounts to using the second equation in Eq.(89) to define the critical temperature. The Jensen-inequality based derivation that shows how this equation limits N is similar to that which leads to Eq.(94). Ultimately, we find whereĒ ≡ D j=1 n j E j /N . Thus the implication is the same as that which follows Eq.(94): Ligandreceptor systems which have a large number of ligands and where optimal binding is physically important, should also have large average binding energies for those ligands.
Eq.(86) gives the thermal condition under which fully optimal ligand-receptor binding can be achieved, but there is not a universal reason for why such a bound configuration would be functionally optimal in all biomolecular systems. Perhaps the microstate where all ligands are in their 8.0 optimal receptor sites is too rigid to be biophysically useful, and thus it is preferred if the system exists in a partially-bound state where most (but not all) of the bound ligands are in their optimal configuration. Thus the value in this formalism might not exist in its ability to predict the temperature at which the system is in a fully bound state, but rather in its power to predict and categorize binding properties above this temperature. For this value, the system-type distinction is the major biophysical utility of this formalism: Ligand-receptor systems are inherently combinatorial, but it seems that if nature were to evolve such systems so as to avoid the kinetic traps of suboptimal bindings, it would engineer ligands to have sufficiently high binding affinities to their optimal-receptors as to easily achieve the combinatorics-limit Eq.(87) at in vivo temperatures.

Discussion
We began this work with the goal of using combinatorics to model how the competition between distinct ligands affects the ligands' equilibrium binding properties to receptors. Before developing a physical model of such binding, we needed to solve a modified version of a well known problem on derangements. After solving this problem and using it to compute a general partition function, we were able to derive implicit formulas for the average number of bound ligands and the average number of optimally bound ligands for each type (Eq.(79)). Narrowing our focus to the case of matched receptor and ligand populations (i.e., r j = n j for all j), we were able to derive the condition under which all ligands are bound to their optimal receptors (Eq.(86)).
The observables defined in Eq.(79) allowed us to compute the equilibrium binding properties of our ligand-receptor system at any temperature (provided we have thermal dependences for γ j and δ j for all j). But the model also provided softer biophysical predictions. Eq.(94) and Eq.(95) can be seen as such predictions, labeled as "soft" because rather than explicitly predicting that some observable has a value, they predict that, in order for a system condition to be satisfied, an observable cannot exceed a certain value. The main implication of these inequalities is that the larger the binding energy advantage for optimal contacts, the more particles the system can have and still be capable of achieving the fully optimally bound state. Thus, systems with more ligands should also have larger optimal binding energies supposing such bindings have functional importance. Note that Eq.(95) was derived for the simple model where free-ligands were taken to be point particles, but an analogous inequality could be obtained for more physically informed values of γ j .
For the task of quantitatively modeling optimal ligand-receptor binding, there are many limitations to the introduced model.
First, we are only partly modeling the combinatorics of our system, a combinatorics that can be characterized by derangements from a pre-defined "correct" or optimal configuration. A more general and flexible model would not a priori define such a configuration and would instead have an interaction matrix defining how various ligand types interact with various receptor types, and would then use this matrix to determine the optimal matchings between ligands and receptors. However, it is not clear how such a more general formalism would be soluble. We started Sec. 3 by representing our problem in terms of a general interaction matrix (in Eq.(13)), but then had to make the simplifying assumption of a pre-defined optimal configuration in order to put the partition function in an analytically tractable form. If we wanted to maintain analytical tractability while pursuing a more general model, we would likely need a similar simplifying assumption.
Second, to derive our various optimal binding conditions (Eq.(57), Eq.(72), and Eq.(82)), we assumed that the system contained the same number of receptor sites as ligands for each type (i.e., n j = r j for all j). This assumption was convenient for the framework of derangements, but does not match a realistic biophysical scenario for which the number of ligands and receptors are not necessarily equal. Thus, a reasonable extension to the Biophysical Implications Sec. 7 would be to find the analytical conditions that define optimal binding for the case of unmatched receptor and ligand populations of each type. However, the general form of the observables in this system (given in Eq.(79)) do not make the assumption of n j = r j for all j, and thus can still be used to numerically compute the number of bound and unbound-particles of various types as a function of temperature.
Third, we did not discuss how the spatial organization of receptors can affect the propensity of ligands to bind to them. Instead we assumed that all receptors were on equal footing in terms of spatial accessibility, and all binding variances could be encoded into the collection of parameters γ = (γ 1 , . . . , γ D ) and δ = (δ 1 , . . . , δ D ). How spatial organization affects binding could likely be better modeled through a matrix-based interaction scheme where spatially occluded receptors have reduced binding affinities, but, as previously mentioned, such an interaction scheme is less analytically tractable than that assumed in the paper. In Sec. 1, we first motivated our initial steps towards the general model presented in Sec. 3 by noting that many ligand-receptor systems exhibit many-to-many interactions in which each type of ligand can bind to many receptors and vice versa. This fact led us to build a model in which even when ligands interacted optimally with a few receptors, these ligands still had the potential to interact with sub-optimal receptors. However, the results of Sec. 6.1 suggest that if optimal ligand-receptor binding (where each ligand binds to its energetically optimal receptor) is functionally important in the system, then that binding has to be highly specific to combat its combinatorial disadvantage. So although our model was motivated by multi-specific ligand binding, the model's results primarily pertain to highly-specific ligand binding and the thermal conditions needed to achieve such binding. Therefore left over from the introduction is the question of how does one quantitatively model the ligand-receptor binding underlying the combinatorial codes that make multi-specific systems so necessary in pharmacology and real biomolecular environments. Building such a model would likely require us to make use of an interaction matrix since it is only through such a matrix that one can encode different sets of ligands interacting optimally with different sets of receptors.
As an ending remark, we note that although Sec. 4.2 explored a limiting case to our more general model, the model introduced there can also stand alone as a probabilistic model of derangements.
For a state space of permutations of a list with repeated elements (e.g., all of the ways to permute 11.0 the characters and spaces in a sentence), we assume that there is a single permutation where the elements are said to be in their "correct order." Let the parameter w i be proportional to the probability that an element of type i is in its correct position. Then, the partition function we previously computed in Eq.(59) becomes the weighted sum of states where G n (defined in Eq. (2) and related to B r,k through G n = B n,n ) is the generalized derangement formula, and n j is n with 1 subtracted from the jth component: n j = (n 1 , . . . , n j − 1, . . . , n D ).
Aside from deriving them directly, we can obtain the expressions in Eq.(97) by translating the Boltzmann factor expressions in Sec. 4.2 into the language of probability weights by taking δ j → w j .
Analogous to Eq.(93), using Jensen's inequality, we find that a necessary, but not sufficient, condition for the average tot to be equal to N (i.e., for all elements to be in their correct order) is This probabilistic model of derangements would be relevant in a computing context where one is interested in the probability of various derangements that differ from a given sequence by a fixed number of elements.

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

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

Author's Statements
• Financial Support: This research received no specific grant from any funding agency, commercial or nonprofit sectors.
• Conflict of Interests Statement: The author has no conflicts of interest to disclose.
• Ethics Statement: This research did not require ethical approval.

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

A Simulations
To simulate our system and create the plots in Fig. 6b, Fig. 7b, Fig. 8b, and Fig. 9 we implemented a Metropolis Hastings algorithm [Kra06]. The code for recreating these figures is linked to in Sec.
12. Here we review the salient parts of the implementation.
For these simulations, we needed to define a microstate, the probability of transitions between microstates, and the types of transitions between microstates.

A.1 Microstate Definition
A microstate of our system was defined by two lists: one representing the collection of unbound particles, and the other representing particles bound to their various binding sites. The particles themselves were denoted by unique strings and came in multiple copies according to the system parameters. For example, a system with D = 3 types of particles with n 1 = 2, n 2 = 3, and n 3 = 1 could have a microstate defined by unbound particles = [A 2 , A 2 , A 3 ] and bound particles = where "−" in the bound list stands for an empty binding site. Since the number of optimally bound particles was an important observable for the system, we also needed to define the optimal binding configuration for the microstates. Such an optimal configuration was chosen at the start of the simulation and was defined as a microstate with no unbound particles and all the bound particles in a particular order. For example, using the previous example, we might define the optimal binding configuration as optimal bound config = [A 1 , A 1 , A 2 , A 2 , A 2 , A 3 ], in which case the number of optimally bound particles of each type in bound particles = [A 1 , −, A 2 , −, A 1 , −] is m 1 = 1, m 2 = 1, and m 3 = 0. The number of bound particles of each type is k 1 = 2, k 2 = 1, and k 3 = 0. We note that the order of the elements in unbound particles is not physically important, but, since the number of optimally bound particles is an important observable, the order of the elements in bound particles is physically important.

A.3 Transition Probability
For these simulations, the temperature-normalized energy of a microstate with k i bound particles of type i and m i optimally bound particles of type i was defined as where k = (k 1 , k 2 , . . . , , k D ), m = (m 1 , m 2 , . . . , m D ), γ i is the binding affinity, and and δ i is the optimal-binding affinity.

A.2 Transition Probability
For transitioning between microstates, we allowed for three different transition types: Particle binding to a site; particle unbinding from a site; permutation of two particles in two different binding sites. Particle binding and unbinding both occur in real physical systems, but permutation of particle positions is unphysical. This latter transition type was included to ensure an efficient-in-time sampling of the state space. For simulations of equilibrium systems it is valid to include physically unrealistic transition types as long as the associated transition probabilities obey detailed balance.
At each time step, we first randomly selected one of the three transition types with equal probability for each type, then randomly selected the final proposed microstate given the initial microstate, and finally computed the probability that said proposal was accepted. By the Metropolis Hastings algorithm [Kra06], the probability that the transition is accepted is given by acceptance prob(init → fin) = min 1, e −β(E fin −Einit) π(fin → init) π(init → fin) , where β is inverse temperature, E init is the energy of the initial microstate state, and E fin is the energy of the final microstate, with energy defined in Eq.(99). The quantity π(init → fin) is the probability of randomly proposing the final microstate state given the initial microstate state and π(fin → init) is defined similarly. The ratio π(fin → init)/π(init → fin) varied for each transition type.

A.3 Transition Types and Examples
Below we give examples of the three types of transitions along with the value of the ratio π(fin → init)/π(init → fin) in each case. In the following, N f and N b represent the number of free particles and the number of bound particles, respectively, before the transition.
• Particle Binding to Site: One particle was randomly chosen from the unbound particles list and placed in a randomly chosen empty site in the bound particles list. π(fin → init)/π(init → Transition weight: π(fin → init)/π(init → fin) = 9/4 • Particle Unbinding from Site: One particle was randomly chosen from the bound particles list and placed in the unbound particles list. π(fin → init)/π(init → fin) = ((N f + 1) −1 × Transition weight: π(fin → init)/π(init → fin) = 3/16 • Particle Permutation: Two particles were randomly selected from the bound particles list, and their positions in the list were switched. π(fin → init)/π(init → fin) = 1. Example: Transition weight: π(fin → init)/π(init → fin) = 1 For impossible transitions (e.g., particle binding when there are no free particles) the probability for accepting the transition was set to zero. At each temperature, the simulation was run for anywhere from 10,000 to 30,000 time steps depending on observed convergence, and the final 2% of the time steps were used to compute ensemble averages of k and m . These simulations were repeated five times, and each point in Fig. 6b, Fig. 7b, Fig. 8b, and Fig. 9 represents the average k and m over these five runs.

B Consistency Checks for B r,k
In this section, we affirm the various consistency checks for B r,k . The first check ensures that B r,k has the proper limiting case when r has components with value 1. The second check ensures that B r−m,k−m has the correct normalization when summed over all possible values of m.
Isolating the summation over the Kronecker delta yields k1 j1=0 · · · k D j D =0 δ(J, j 1 + · · · + j D ) = 1 2πi Thus, we obtain Next, we make the change of variables We then have To simplify this result we need to find an identity for Expanding the generalized Laguerre polynomial L (α) q (x) according to its definition, we obtain More identity wrangling gives us and thus Eq.(117) becomes which, by the definition of the Laguerre polynomial yields Also, from the definition of the Laguerre polynomial, we can show Therefore, with Eq.(120) and Eq.(121), we have k q=0 k + α k − q (−1) q L (α) q (x) = lim Inserting this result into Eq.(115) gives us which, with the Gamma function definition, yields I r,k = (r 1 + · · · + r D )! (r 1 − k 1 + · · · + r D − k D )!k 1 ! · · · k D ! .
Why does this result make sense? From one perspective, the quantity r1 m1 · · · r D m D B r−m,k−m is the number of ways to choose m j out of r j positions to contain their correct elements while the remaining k j −m j elements are completely deranged with respect to their r j −m j remaining correct positions, for j = 1, . . . , D. If we sum over all possible values of m j we should obtain the number of ways to arrange (i.e., not only derange) k j ≤ r j objects for j = 1, . . . , D across a total of r 1 + . . . + r D lattice sites.
On the other hand, if we are trying to calculate the number of ways to arrange k j objects of type j for j = 1, . . . , D across r 1 +· · ·+r D lattice sites, we can use the language of multinomials. Say that the objects represent "filled" sites on the lattice and the spaces between objects are "empty" sites. Then there are k 1 + · · · + k D filled sites (of which k i are identical for each i) and r 1 − k 1 + · · · + r D − k D empty sites. Finding the number of ways to arrange the objects amongst the r 1 + · · · + r D sites is equivalent to finding the total number of ways to order this collection of filled and empty sites while correcting for equivalent orderings due to reordering the positions of the same type of site.
Including both filled and empty sites there is a total of r 1 + · · · + r D sites amongst which we have r 1 − k 1 + · · · + r D − k D "copies" of empty sites, k 1 copies of filled sites of type 1, k 2 copies of filled sites of type 2, . . . and k D copies of filled sites of type D. Counting the number of ways to order the r 1 + · · · + r D sites and correcting for the equivalent reorderings arising from the multiple copies of various types of sites leads to the multionomial (r 1 + · · · + r D )! (r 1 − k 1 + · · · + r D − k D )!k 1 ! · · · k D ! , which we have shown is equivalent to the result written in terms of a summation over B r−m,k−m . C.0

C Deriving General Partition Function
In this section, we derive Eq.(24) from Eq.(22). First we note from the definition of B r,k in Eq. (6) that B r−m,k−m = 1 ( j α j )!
where we changed variables from m j to q j = k j − m j in the final line. Using the following Laguerre polynomial identity (derived in Eq.(119)) k q=0 k + α k − q U q L (α) q (X) = (1 + U ) k L (α) k XU 1 + U , we then have x αj (δ j − 1) kj L (αj ) kj For the summation over k, we use the contour integral identity for inverse factorial

D Laguerre Polynomial Argument Summation Identity Eq.(133)
where Γ is a closed contour about the origin, to obtain k C r,k n + r − n n − j where we changed variables to j = n − in the fourth line and used the definition Eq.(5) in the final line.
Next, we seek to eliminate thez andx dependence from these observables. From Eq.(137), Eq.(138), and Eq.(139), we can showz where we used k = D j=1 k j , and the total particle number N is equal to the number of particle types D when there is one copy per particle. With Eq.(138) and Eq.(139), we can also show which is our first equilibrium condition.
which is our second equilibrium condition.

F Derivation of Eq.(82)
In this section, we derive Eq.(82), the general thermal condition for the fully optimally bound state.
This state is defined as where N ≡ D j=1 n j = D j=1 r j . To derive the desired thermal condition, we will need some of our previous results. In particular, we need the average number of optimally bound ligands of type j m j = n j δ j δ j − 1 and the equations relatingz andx to our observables: From Eq.(147), we can infer that for the fully optimally bound state we have where we are explicitly referencing the O(D −1 ) terms in order to satisfy where we subsumed terms of order O(D −1 δ −1 ) into O(δ −1 ).
For the fully optimally bound configuration, ligands need to greatly favor their optimal receptors in order to bind only to such receptors. Thus for the fully optimally bound state, we can assume that the optimal-binding affinity for each particle-type is much greater than 1: δ j 1. From this assumption we can takeφ j defined asφ to satisfy |φ j | 1. We will check this latter assumption at the end of the calculation, but moving forward with it, we find With k = N − 1 + O(δ −1 ) andz = N − k , we findz = 1 + O(δ −1 ), thus giving us the final thermal