Introduction

For an aggregation to take place between allergen, DARPin, and IgE molecules, these molecules must first go through a complexation phase. Complexation consists in non-covalent binding of IgE to allergen and/or DARPin molecules. It is a necessary step which thus controls aggregation. A complexation model was thus built to maximize the chances of aggregation. As shown in Figure 2, a too low concentration of cells (hence of allergen and DARPin molecules) will cause saturation by IgE and then prevent any aggregation; inversely, if too many cells are present the aggregation will not take place either. We thus developed a model to determine the concentrations of allergen and DARPin molecules which maximize the pool of species capable of interacting to form an aggregation (Figure 2). These optimal concentrations were calculated for a wide range of allergens (defined by their dissociation constant, KD), to ensure our detection system will be as universal as possible to cover the complete range of allergens.

Figure 2: Two extreme cases of IgE concentration preventing aggregation. As an example for a constant number of IgE, here 16, the presence of aggregation or not is directed by the number of strains present in the system. If there is not enough bacteria in comparison to the number of IgE, all sites will be saturated, and there will not be aggregation (A). The maximum aggregation corresponds to an optimal repartition of IgE, where IgE can bind both strain A and D (B). If the concentration of bacteria is too high, IgE will bind, but there will not be enough IgE to make a consistent aggregation.

The reaction of complexation

The complexation phenomenon, prior to aggregation, takes place between allergen (A), DARPin (D), and IgE molecules (Y). As the constant part of the IgE as a high affinity for the DARPin and its variable part a high affinity for the allergen, these molecules can bind and form two types of complexes: allergen linked to an IgE (AY) and DARPin linked to an IgE (DY). One should note that complexation is a reversible mechanism. The corresponding reactions are schematized in Figure 3.

Figure 3: Reactions of complexation forming the AY and DY complexes. Strain A reacts with Y (i.e, IgE) following the rate constant \(r_{1}\). Strain D reacts with Y following the rate constant \(r_{2}\).

The rates of these association/dissociation reactions can be modeled using the reversible law of mass action with two rate constants:

\( r_{1}=k_{on_{1}}[D][Y]-k_{off_{1}}[DY]\qquad(1)\)

\( r_{2}=k_{on_{2}}[A][Y]-k_{off_{2}}[AY]\qquad(2)\)

The \(K_{D}\) of allergen and DARPin were more documented than the corresponding \(k_{off}\) (Table 1), so we introduced this parameter in the above equations, which yield:

\( r_{1}=k_{on_{1}}([D][Y]-[DY]K_{D_{DY}})\qquad(3)\)

\( r_{2}=k_{on_{2}}([A][Y]-[AY]K_{D_{AY}})\qquad(4)\)

With \( K_{D}=\frac{k_{off}}{k_{on}}\qquad(5)\)

From these equations we obtained the ODE system describing the evolution of the concentrations over time:

\( \frac{d[D]}{dt}=-r_{1}\qquad(6)\)

\( \frac{d[A]}{dt}=-r_{2}\qquad(7)\)

\( \frac{d[Y]}{dt}=-r_{1}-r_{2}\qquad(8)\)

\( \frac{d[DY]}{dt}=r_{1}\qquad(9)\)

\( \frac{d[AY]}{dt}=r_{2}\qquad(10)\)

We focused on the steady-state of this system, when an equilibrium between each species is finally reached. Note that the steady-state depends only on the \(K_{D}\) of each reaction and is independent of the rate constants, so we could fix an arbitrary value for the \(k_{on}\) of both reactions.
These ODE were implemented in the open-source software COPASI (v4.36) and solved by its integrated ODE solver. All simulations were carried out via the CoRC (v0.11.0) package of R (v4.0.2). You can find our programs used for the aggregation model in the Gitlab folder: Complexation model. All parameters are listed in Table 1.

Conversion between molecules and cell concentration

As we want to determine the optimal concentration of each strains, we converted the concentrations of molecules into the corresponding concentration of cells (expressed as optical density, more appropriate and easier to compare with wet lab results), using the following formula:

\( OD_{600_{nm}}=\frac{N_{a}}{X}*\frac{[A]+[AY]}{n}\qquad(11)\)

where \(N_{a}\) is the number of Avogadro, \(X\) is the concentration of E. coli in cell/mL corresponding to an \(OD_{600_{nm}}\) of 1, and \(n\) is the total number of A-molecules displayed on the surface of each cell. Importantly, this parameter depends on the expression level, which can be useful to further optimization. These parameters are listed in Table 1.

Parameters

As we wanted to cover a large range of allergens, we chose Allergen - IgE \(K_{D}\) ranging from \(10^{-11}\) to \(10^{-7}\,M\). In that way, we have chosen allergens showing representative \(K_{D}\) values of this range. Values of these \(K_{D_{AY}}\) are listed in Table 1.

Table 1: Parameters of the complexation model. This table includes the parameters used for the conversion in \(OD_{600_{nm}}\) as well as those of the reactions of complexation.

Initial conditions

The initial concentrations of A and D molecules were chosen so that the corresponding \(OD_{600_{nm}}\) equals 0.1, which is high enough for easy detection on our spectrophotometer but supposedly not too high to generate non-specific aggregates.
We chose the concentration of IgE as the minimal concentration at which the test results must be positive. We defined this value as the cut-off value of Immunocap tests (i.e. the minimum value to diagnose a patient as allergic) (Dang et al., 2012). As the cut-off value corresponding to an allergic patient is greatly different depending on the nature of the IgE, we chose to calculate the value for anti-Ara h 2 IgE only, the IgE used for our aggregation assays. The literature indicates a cutoff value of < 0.15 kU/L for patients with a peanut allergy. According to Merck manuals, it equals 0.4 µg/L, which equals \(2.1*10^{-12}\,M\). We considered that the antibodies were diluted in the reaction mix, with a dilution ratio of 1:10 (as could be done in a real situation), so their initial concentration was set as \(2.1*10^{-13}\,M\).
Initial concentrations of all species are listed in Table 2.
If we consider there are no complexes at t=0, the initial concentrations of AY- and DY-complexes are set to zero.

Table 2: Initial concentrations of all species. The initial concentrations of complexes are set to zero.

Potential of aggregation

Once the system has reached a steady-state (i.e. the concentration of free and complexed molecules are stable over time), the couples of particles capable of interacting to form an aggregation are A+DY and D+AY, while the other couples of particles (A+D and AY+DY) cannot interact. We thus defined a metric to quantify the potential of aggregation as a function of the concentration of each particle.
We take as example A and DY to explain our approach. First, \(min(A,DY)\), i.e. the minimum between the number of A and DY particles, represents the maximal possible number of interactions that can initiate aggregation. For example, if three A-molecules are put in the presence of five DY-complexes, 3 interactions capable of forming an aggregation can occur at most.
Then, we introduced the terms \(\frac{A}{A+AY}\) and \(\frac{DY}{D+DY}\), which are respectively the fraction of free A-molecules among the total pool of allergen molecules, and the fraction of DY-complexes among the total pool of DARPin molecules. Assuming homogeneous distribution of the free and complex molecules on the surface of all cells, these terms represent the proportion of allergen and DARPin molecules on the surface of each strain that may initiate aggregation. A high proportion of molecules in the correct complexation state (A and DY) is expected to increase the chance for the cells to bind, and finally aggregate. In contrast, a low proportion of these molecules will reduce the strength of interaction between the cells, and reduce the chance of aggregation.
Therefore, we defined the aggregation potential of the couple A+DY as the product between the maximal number of possible interactions for this couple and the chance for these interactions to occur when two cells collide, and thus initiate complexation:

\( Aggregation\,potential(A+DY)=min(A,DY)*\frac{A}{A+AY}*\frac{DY}{D+DY}\qquad(12)\)

The same reasoning applies for AY+D, hence the aggregation potential of the complete system is defined as:

\( Aggregation\,potential=min(A,DY)*\frac{A}{A+AY}*\frac{DY}{D+DY}+min(AY,D)*\frac{AY}{A+AY}*\frac{D}{D+DY}\qquad(13)\)

This metric has the following properties:

• - If there is an excess of cells compared to IgE, the fractions \(\frac{AY}{A+AY}\) and \(\frac{DY}{D+DY}\) tend to 0, so the aggregation potential also tends to 0: there is no interaction. Note that the same situation is encountered when no (or very few) IgE are present.
• - If there are much less cells than IgE (all sites saturated by IgE), the fractions \(\frac{A}{A+AY}\) and \(\frac{D}{D+DY}\) tend to 0, so the aggregation potential also tends to 0: there is no interaction.

These two cases represent a saturation of cells (case 1) and of IgE (case 2). In these situations, the aggregation potential falls to 0, as expected.

Simulations for optimal experimental design

We then used the model to determine the aggregation potential in silico for different conditions that could be encountered in the DAISY laboratory. Our approach was illustrated with one specific allergen: Ara h 2. It was then extended to a wide range of allergens (with a \(K_{D_{AY}}\) range of \(10^{-11}\) to \(10^{-7}\,M\)), to determine the optimal concentrations of A and D molecules for as many allergens as possible.

A first simulation was made: calculated aggregation potential as a function of concentration of A-molecules. The plotted results are shown in Figure 4.

Figure 4: Aggregation potential calculated for different initial concentrations of A-molecules. The A molecule concentrations range from \(10^{-14}\) to \(10^{-7}\,M\) (the other initial concentrations being as in Table 2). The simulation was done with 50 steps.

The aggregation potential reached a maximum at \(10^{-11}\,M\), corresponding to the optimal initial concentration of A for this particular case. When the concentration of A-molecules is too low the number of A and AY approach zero, so are \(min(A,DY)\) and \(min(AY,D)\) and then the aggregation potential also tends to zero. Similarly, when their concentration is too high the IgE are too diluted, so the fractions \(\frac{AY}{A+AY}\) and \(\frac{DY}{D+DY}\) and then the aggregation potential tend to zero.

The same simulation was done but this time varying both A and D concentrations, and calculating the corresponding aggregation potential. Several allergens were tested, differing by their \(K_{D}\). The results plotted on contour graphs are shown in Figure 5.

Figure 5: Evolution of the aggregation potential for different initial concentrations of A and D free molecules, for different allergens: Ara h 2 (A), Fel d 1 (B), Ana o 3 (C), Der p 2 F (D), and Der p 2 H10:H12 (E). The A and D molecule concentrations range from \(10^{-14}\) M to \(10^{-7}\) M. The simulations were done with 50 steps. Note that the color skin for the aggregation potential changes from one panel to another.

The maximal value of the aggregation potential was higher for the allergens with low \(K_{D}\). The optimal concentrations varied greatly depending on the allergen. More precisely, the profile of concentrations changed between \(K_{D}\) inferior and superior to \(K_{D_{DY}}\) (equal to \(6.29*10^{-9}\,M\)): a vertical profile for \(K_{D_{AY}} < K_{D_{DY}}\), an horizontal one for \(K_{D_{AY}} > K_{D_{DY}}\). We then decided to calculate the aggregation potential for a range of \(K_{D}\) comprising those of our allergens: from \(10^{-11}\) to \(10^{-7}\,M\). This range has a magnitude 100 times larger and 100 times smaller than the \(K_{D_{DY}}\), so it was representative of the evolution of the aggregation potential. These results were finally compiled into one plot of optimal concentrations versus Allergen - IgE \(K_{D}\), as shown on Figure 6.

Figure 6: Ranges of optimal initial concentrations of A- and D-molecules depending on Allergen - IgE \(K_{D}\). The ranges correspond to the concentrations for which the aggregation potential is at more than 80% of its maximal value.

Ranges of optimal concentrations appeared to change depending on the Allergen - IgE \(K_{D}\). For \(K_{D}\) below \(10^{-9}\,M\), the range of optimal A-molecules was very narrow (more or less than \(10^{-11}\,M\)) while that of D-molecules was quite large (from \(10^{-13}\) to \(10^{-10}\,M\)). For \(K_{D}\) over \(10^{-8}\,M\) the trend reversed: the range of optimal A-molecules became very large (\(10^{-13}\) to \(10^{-8}\,M\)) while that of D-molecules narrowed (around \(10^{-8}\,M\)). For 100,000 protein sites at the surface of bacteria, it corresponds to \(OD_{600_{nm}}\) from \(6*10^{-2}\) (for \(10^{-8}\,M\)) to \(6*10^{-6}\) (for \(10^{-13}\,M\)). These data will be useful in the wetlab as they give ideal initial concentrations of each strain to maximize aggregation, and this for a large panel of allergens.

Conclusion

The first round of our modeling effort has shown useful results for our aggregation assays on the wetlab. This model has shown that aggregation is a quite robust phenomenon. Indeed, optimal ranges of concentrations allowing aggregation have proven to be large enough for an experimental use, but we probably did not use the accurate concentrations of strains A and D during our experiments (they were too concentrated), explaining our difficulties to aggregate the strains. Results of this modeling effort with the complexation model have been a key tool to improve our aggregation assays, so they will be integrated on our DBTL cycle to improve future tries of aggregation (see Engineering success). The complexation model has brought us a deeper understanding of the impact of concentrations and the expression of DARPin or allergens on the bacterial membrane. With this model, we have been able to determine the optimal concentrations for allergen with any affinity (\(K_{D}\)), as well as a range of concentrations for which an unknown allergen should give rise to aggregation. This represents an important step towards our versatile detection system able to diagnose any allergy.

Introduction

With the complexation model, we highlighted the deep impact of each element concentration on the strength of interactions, and thus on the potential of aggregation. However, aggregation is a dynamic process. To maximize the speed of aggregation (and thus indirectly the throughput-compliancy of our method), we developed a dynamic model of cellular aggregation. This dynamic model will help us to consider the molecular scale to integrate other experimental set-up for our aggregation assays. However, a very complex model is required to understand the time needed and the final size of aggregates, so this model will be built with an iterative process between the dry lab and the wet lab. This will help us to design a relevant model, by starting with simple equations and strong hypotheses, to give a first direction of the set-up for the wetlab. Data produced by aggregation assays will help us to improve the model. We describe in details the construction of the model and the first loop of our iterative improvement cycle in this section.

Aggregation dynamics model

To simulate aggregation, we developed a population balance equation (PBE) model. Basically, this modeling approach considers that four processes define the behavior of any aggregate (Jeldres et al., 2017):

• - A birth by aggregation of two smaller particles to make the particle considered (\(B_{a}\))
• - A death by aggregation with another particle to obtain a bigger aggregate (\(D_{a}\))
• - A birth by break of another bigger aggregate (\(B_{b}\))
• - A death by break of the aggregate considered into two smaller particles (\(D_{b}\))

These processes are shown in Figure 7.

Figure 7: representation of each term \(B_{a}\), \(D_{a}\), \(B_{b}\), \(D_{b}\) for the given particle of size N=2. \(B_{a}\), \(D_{a}\), \(B_{b}\), \(D_{b}\) are defined for the number of cells on particles N=2.

Here, we consider two different particles: the cell displaying an allergen called cell A, and the cell displaying DARPin named cell D, considering that the complexation phenomenon with IgE is already done with the maximum potential of interaction. We thus designed a bivariate PBE (Fox, 2006). The two variates used here are \(n_{A}\) and \(n_{D}\) representing both the number of particles A and D on an aggregate. The dynamic of any aggregate can be represented by the following equations:

\( \frac{\partial f(n_{A}, n_{D})}{\partial t}=B_{a}-D_{a}+B_{b}-B_{b}\qquad(14)\)

The aggregation phenomenon can be separated physically into two phases. As it happens, during the beginning of the aggregation phenomenon, there is aggregation with single cells, then aggregation between aggregates. Over time, the size of aggregates is increasing. Then during the second phase, aggregates break when they become too big to be stable or the velocity gradient is too high (Wilson et al., 2019).
Therefore, we will build a first version of the PBE model considering only the first phase corresponding to the aggregation phenomenon and use soft agitation to avoid disruption of our aggregates so that we could neglect breaking effects in our PBE.
Based on the literature, we expect to detect aggregates before reaching a critical size of the order of a millimeter where the break phenomenon could not be avoided (Kylilis et al., 2019). This hypothesis allows us to build a first PBE model considering only aggregation terms, and neglecting breaking effects \(B_{b}\) and \(D_{b}\) to study the first phase of aggregation, yielding the following system of equations:

\( \frac{\partial f(n_{A}, n_{D})}{\partial t}=B_{a}-D_{a}\qquad(15)\)

The \(B_{a}\) term corresponds to the birth of an aggregate resulting in clumps of two smaller aggregates. This term can be modeled using the following equation that we explain in detail below:

\( B_{a}=\frac{1}{2}\int_{1}^{n_{D}-1}\int_{1}^{n_{A}-1}f(n'_{A},n'_{D})*f(n_{A}-n'_{A},n_{D}-n'_{D})*\beta(n'_{A},n_{A}-n'_{A},n'_{D},n_{D}-n'_{D})*\Gamma(n'_{A},n_{A}-n'_{A},n'_{D},n_{D}-n'_{D})*dn'_{A}dn'_{D}\qquad(16)\)

Here we have:

• - \(n'_{A}\): the number of strains A in aggregates of any size
• - \(n'_{D}\): the number of strains D in aggregates of any size

The formation of aggregates of a given size is proportional to the number of smaller aggregates multiplied by a collision frequency and a collision efficiency. To consider all possibilities of clumping aggregates, this product of terms is integrated over all aggregates with a size smaller than (\(n_{A}\), \(n_{D}\)). The half ponderation results from the fact that the integral term naturally considers each collision event twice (by symmetry).

This equation considers a clump of two complementary aggregates represented by the \(f(n'_{A},n'_{D})\) and \(f(n_{A}-n'_{A},n_{D}-n'_{D})\). The product of these proportions gives the proportion of the new aggregate (Fox, 2006).

To direct aggregation, two laws \(\beta\) and \(\Gamma\) have to be formulated. The first law defines the rate of encounter as the number of collisions per unit time. In general, this collision frequency is scaled by the velocity gradient (Meyer et al., 2011). Unfortunately, we do not have yet the appropriate aggregation experiment to choose the one physically relevant for our model, so we chose a standard formulation, proportional to the mean of velocity gradient:

\( \beta(n'_{A},n_{A},n'_{D},n_{D}=K_{\beta}*G)\qquad(17)\)

with:

• - \(K_{\beta}\): probability of efficient encounter
• - \(G\): mean gradient of hydrodynamic velocity of the system

The second law, \(\Gamma\), defines the efficiency of the collision (it is actually the formation of new aggregates after collision). Here we had to design our own proposal for the aggregation efficiency between different strains (see Demonstration of the \(\Gamma\) law), considering the number of DARPin \(U_{SD}\) and allergens \(U_{SA}\) on the surface. The product of these two terms gave the factor of collision efficiency \(K_{\Gamma}\):

\( K_{\Gamma}=U_{SA}*U_{SD}\qquad(18)\)

So the \(\Gamma\) law becomes:

\( \Gamma(n'_{A},n_{A},n'_{D},n_{D})=K_{\Gamma}\frac{n_{A}n'_{D}+n_{D}n'_{A}}{(n_{A}+n'_{A})(n_{D}+n'_{D})}\qquad(19)\)

Then, we had to look to the \(D_{a}\) term, the death of aggregates of the considered size, resulting from a clump of two aggregates to form a bigger one:

\( D_{a}=f(n_{A},n_{D})\int_{1}^{\infty}\int_{1}^{\infty}\beta(n'_{A},n_{A},n'_{D},n_{D})*\Gamma(n'_{A},n_{A},n'_{D},n_{D})*f(n'_{A},n'_{D})*dn'_{A}dn'_{D}\qquad(20)\)

The death term describes the aggregation of one particle of considered size and any other particles using the same \(\beta\) and \(\Gamma\) laws (Fox, 2006).

Finally, the aggregation problem is formulated as:

\( \frac{\partial f(n_{A}, n_{D})}{\partial t}=\frac{1}{2}\int_{1}^{n_{D}-1}\int_{1}^{n_{A}-1}f(n'_{A},n'_{D})*f(n_{A}-n'_{A},n_{D}-n'_{D})*\beta(n'_{A},n_{A}-n'_{A},n'_{D},n_{D}-n'_{D})*\Gamma(n'_{A},n_{A}-n'_{A},n'_{D},n_{D}-n'_{D})*dn'_{A}dn'_{D}\)
\( -f(n_{A},n_{D})\int_{1}^{\infty}\int_{1}^{\infty}\beta(n'_{A},n_{A},n'_{D},n_{D})*\Gamma(n'_{A},n_{A},n'_{D},n_{D})*f(n'_{A},n'_{D})*dn'_{A}dn'_{D}\qquad(21)\)

With \(\beta\) and \(\Gamma\) defined in equations (17) and (18) respectively.

Implementation

The PBE equation was implemented in Matlab and solved using the DQMOM method (Marchisio et al., 2005; Fox, 2006) (see Presentation of the DQMOM). Briefly, the general principle of this method consists in choosing some representative values of distribution function \(f(n_{A},n_{B})\) by choosing two abscissas representing both \(n_{A}\) and \(n_{D}\) and ponderate them with weights. This solving method then simulates the evolution of moments, a magnitude integrating abscissas and weights over the time. This approach allows to simplify the model by solving only a few equations are solved in place to the full distribution \(f(n_{A},n_{B})\). Consequently, the DQMOM method is computationally swiftest. The main script SimulationX initializes parameters, and solves the differential equation with a differential equation solver. This solver calls the function DQMOM_Solver, a function able to solve with the DQMOM method the PBE model written on PBE_Model. You can find our programs used for the aggregation model in the Gitlab folder: Aggregation model.

Thanks to the complexation model described above, we know the optimal number of DARPin and Ara h 2 allergen needed on the surface and the optimal concentration of strain A and D. The last critical set-up parameter is the agitation represented by the velocity gradient \(G\) on the aggregation dynamics model. So we conducted one simulation experiment to show the relevancy of the model to predict the dynamic of aggregation and to test the impact of the velocity gradient for a certain timescale, here 1 hour. Without aggregation assays on the wet lab, the \(K_{β}\) cannot be evaluated, so an arbitrary value was used, for that 0.01 was chosen. Three values of velocity gradient (5, 10 and 20 \(s^{-1}\)) were used to run simulations to have an idea of the impact of the velocity gradient. This scale of velocity gradient corresponds to a usual soft agitation on a 2 mL tube (Liné, personal communication, August 1, 2022).

Parameters values used in the first version of the model are given in Table 3.

Table 3: parameters of the aggregation model.

Results and interpretation

First, for the velocity gradient \(G=10\;s^{-1}\) the relevancy of our model was tested by simulating the total amount of each strain A and D (aggregated or not), hence verifying the mass conservation of cells in the system (Figure 8).

Figure 8: Validation of the relevancy of the resolution by the DQMOM method. The proportion of each strain A and D on the system was followed during the simulation. This validates the parameters used for the DQMOM method.

As shown in Figure 8, the amount of cells A and D are conserved over time. This validates the DQMOM method and its implementation (see Presentation of the DQMOM), which can then be used to test the impact of some parameters, here the agitation, on the aggregation dynamics (Figure 9).

Figure 9: Aggregation dynamics evolution as function of the velocity gradient. The mean number of aggregates was calculated with the DQMOM method to compare the impact of the velocity gradient (i.e. agitation).

We thus used the model to run different simulations by varying the agitation, represented here by the velocity gradient \(G\). Results confirm our intuition that when the velocity gradient increases, it improves our aggregation dynamics, and raises the final size of aggregates as a consequence of a higher collision frequency. According to this first version of the model, the dynamics of aggregation is incredibly fast and strong: in less than 100 seconds, aggregates reach a mean number of 1,000 cells. After 500 seconds, the size of aggregates tends to be constant, there is a limitation on the aggregation phenomenon, but it continues to increase since the breaking effects are not taken into account to counterbalance the aggregation. Implementing this breaking phenomenon in the next version of the model is expected to slow down the agregation kinetics. For aggregation assays on wet lab, using a stronger agitation should help us to make aggregates bigger by increasing the chance of particles to encounter. These assays will have to be carried out, and we will use the corresponding data us to refine our aggregation terms on our PBE and estimate parameters values.

Conclusion

With our strategy, we were able to validate the structure of the agregation model and of the simulation algorithm, and to simulate in silico a first aggregation experiment. According to the first version of this model, agitation is a key parameter for optimizing agregation.
The next steps will consist in conducting wet lab experiments to monitor the size of aggregates along time. Following the DBTL cycle, these data are now required to refine and improve the model, in particular by:

• - Validating the relevancy of the aggregation law implemented here, or choose an alternative law which matches better the physics of aggregation which is specific to the type of particules considered.
• - Determine the value of each parameters (i.e. calibrating the model) by fitting experimental kinetics of aggregation measured in the wet lab experiments.
• - Implementing the breaking phenomenon, which will allow us to improve the accuracy of the predicted aggregation dynamics.