# Simulating Mono-and Multi-Protein Phosphorylation within Nanoclusters

Olivier Destaing, Bertrand Fourcade

► To cite this version:

Olivier Destaing, Bertrand Fourcade. Simulating Mono-and Multi-Protein Phosphorylation within Nanoclusters. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2025, 10.1103/bb12-kz72. hal-05356897

HAL Id: hal-05356897

https://hal.science/hal-05356897v1

Submitted on 10 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Simulating Mono- and Multi-Protein Phosphorylation within Nanoclusters

Olivier Destaing

Institute for Advanced Biosciences,

Centre de Recherche Université Grenoble Alpes,

Inserm U1209, CNRS UMR 5309, France.

Bertrand Fourcade

Université Grenoble Alpes and CNRS, LIPHY, F-38000 Grenoble, France

(*Bertrand.Fourcade@univ-grenoble-alpes.fr)

(Dated: November 10, 2025)

Protein nanoclustering is a characteristic feature of their activated state and is essential for forming numerous subcellular structures. The formation of these nanoclusters is highly dependent on a series of post-translational modifications, such as mono- and multi-phosphorylation and dephosphorylation of residues. We theoretically simulate how a protein can be either mono- or multi-phosphorylated on several residues in functional nanoclusters, depending on effective biophysical parameters (diffusion, dwell time, etc.). Moving beyond a binary view of phosphorylation, this approach highlights the interplay between mono- and multi-phosphorylation, the cooperative effects generally associated with multi-phosphorylation networks, and stresses the role of phosphatases in transforming graded phosphorylation signals into almost switch-like responses. The results are discussed in light of experiments that probe the distribution of phospho-residues.

# I. INTRODUCTION

Living cells process information from their extracellular environment through biochemical signaling networks, where protein phosphorylation is one of the mechanisms responsible for altering their biological activity, location, or docking with other proteins [1]. These regulatory mechanisms are orchesterated by a large set of kinases and phosphatases, where a protein can be mono- or multi-phosphorylated on different residues [2, 3]. Protein phosphorylation is a local and contextual phenomenon, meaning the same system of molecular components can display different properties in different situations.

An emerging concept in many areas of biology, including cell and developmental biology or immunology, is the existence of finite-lifetime platforms at sub-resolution length scales where proteins are phosphorylated [4–10]. This non-homogeneous spatial organization aligns with reports showing 40-250 nm nanoclusters of different adhosome proteins throughout adhesions, which organize these protein covalent modifications for signal and information transmission [11, 12].

Following classical works that established the first-passage framework [13–15], this paper reformulates these concepts within reaction-diffusion models for nanoclusters to show how essential properties can be combined [9, 16]. The first, called pleiotropy, refers to the possibility for the same signaling element—a kinase, phosphatase, or GTP-ase, for example—to induce different responses. The second is the possibility of crosstalk between signaling pathways, where one signaling input is amplified or switched off depending on a second one. The third is the specific decision-making that explains how a single signaling event can generate different cellular responses. Although discussed in the biological literature since the 1980s and in physics more recently [17], the decision-making process between antagonistic pathways has mainly been studied at the molecular level in a binary representation

where a protein undergoes post-translational modifications [18]. While important, these works provide only a limited understanding of experiments where, for example, optogenetic control of signal processing modifies and selects cellular responses [19–23]. Model systems that combine molecular data with biophysical constraints could therefore be useful in a context where changing the dynamics of a single molecule elicits drastic changes, as demonstrated by spatio-temporally resolved phospho-proteomics [24, 25].

As applied here, bidimensional ligand-receptor networks at the front-end of cell signaling system are complex signal-units able to select among different pathways for different environment [26, 27]. Since phosphorylation by kinases is a key regulatory layer for protein-protein interaction networks, these experiments and these conceptual approaches prompt the inclusion of 2D diffusion effects. By including diffusion into effective nanoscopic phosphorylation rates, the aim here is to show, from the perspective of mono-and-multiphosphorylation, how the following questions can be answered: How is a substrate selected? What is the probability for a protein to be phosphorylated? What is the phosphorylation mode, and how does its likelihood change with varying parameters? Is it an all-or-nothing process, an almost switch-like process, or a graded rheostat-type process? Answering these questions is a first step toward a better understanding of a system capable of different molecular self-organization using the same set of proteins [28, 29]. We will answer this question here by calculating exactly the 2D diffusive flux of a molecule with different boundary conditions, as this problem has already been addressed in Ref. [30].

Biophysical models for mono- and multi-phosphorylation in signaling to model information processing have been the subject of previous works based on different versions of mass-action kinetics [31–34]. They are generally based on the McKeithan model [35], which generalizes the Hopfield-Ninio paradigm [36, 37]. In these models, an enzyme-bound substrate must complete a series of phosphorylation events for
a cellular response to occur, see [38] and references therein. Here, the spirit is the same, but the molecular kinetic constants are replaced by effective parameters that allow the system to probe biophysical constraints, such as the size of clusters, the dwell time of proteins specifically anchored on the plasma membrane, or their diffusion, which can vary from one zone to another. Our approach differs from that in other references, such as [39, 40], which assume an escape probability for diffusion to modify the kinetic rate constants in 3D. Since random walks are recurrent in 2D [15], we take a different approach and we solve the full diffusion time-dependent diffusion problem from an analytical and stochastic point of view. Crosstalks between different pathways can then be naturally integrated into this approach where one signaling pathway reads local, averaged biophysical properties modified by another. This spatiotemporal approach also takes into account diffusive and rebinding effects, which are known to induce significant changes in molecular complexes produced in response to the phosphorylation of substrates, leading to changes in the signaling network [41–44] or to the induction of biomolecular condensates [45–48].

This paper focuses on the calculation of the probability for a protein to be mono- or multi-phosphorylated. Multi-phosphorylation is crucial since it can involve a large number of states whose stability may depend on this number [49, 50]. Thanks to its combinatorial properties, the present paper shows that the process of multi-phosphorylation combined with its reverse process leads to very different responses depending on the averaged local environment, reflecting the presence of different signaling pathways. Multi-phosphorylation also allows for cooperativity. In cell signaling, cooperativity is often defined as the result of an allosteric phenomenon, where a chemical change in one of a protein's residues (in this case, phosphorylation) affects the affinity of another residue for another reaction. However, multiple binding can also be induced by multi-phosphorylation-sustained cooperativity, a process often referred to as avidity [51], which we know can lead to the phenomenon of superselectivity due to entropic factors [52–54]. Because of this avidity effect, bivalent induced phospho-binding with proper spacing, even contributed by two weak bonds, can produce stable protein complexes and may be important for understanding the specificity of some signaling phosphatases [55, 56]. We show here that phosphatases, when integrated into the model, have strong effects on the calculation of these probabilities, with potential switch-like effects.

In short, for a relatively small number of residues, the model presented below demonstrates that multi-phosphorylation leads to new phenomena including, for example, cooperativity, changes in critical time scales, and graded versus almost switch-like phosphorylation modes. Steep variations and abrupt changes in phosphorylation likelihoods have the potential to have consequences on the biological side, and this work highlights the importance of having multiple interaction sites for phosphorylated substrates, as described by the modular logic for some kinases and phosphatases [55].

We proceed both analytically, within a type of effective

medium approximation, and numerically. Numerical simulations using a spatial Gillespie's algorithm [57, 58] are useful because they provide access to the stochastic regime for systems with very few molecules and to quantities that are difficult to obtain analytically. The method is standard and is briefly summarized in Appendix A. This paper is structured as follows. The first sections introduce the notations and specify the various models. The case of mono-phosphorylation will then be discussed. Multi-phosphorylation and, finally, the effect of phosphatases will be discussed last. The conclusion summarizes these elements and briefly addresses the potential biological applications, while technical details can be found in the appendices.

## II. THE MODELS

In parallel with Single Particle Tracking experiments that probe the journey of a fluorescent protein within integrin-based structures (see [59, 60], for example), we model a process where a protein is deposited at time $t = 0$ anywhere within an activation disk of size $R_a$. Efficient binding to a partner requires mono- or multi-phosphorylation of key residues. The probability of success of this process and its average time are assumed to be measurable. To simulate this process, we attach a virtual counter to this protein and stop the process when a given number $n$ of residues is phosphorylated, with $n = 1$ for mono-phosphorylation and $n > 1$ for multi-phosphorylation. Having a given number of phospho-residues is a specific outcome out of several alternative possibilities, since the protein can desorb from the membrane and return to the cytosol or be absorbed on specific boundaries. This is an example of splitting probability where the process can split into two different outcomes [15].

The physical framework presented here treats free diffusion along the membrane and efficient phosphorylation as two generic processes. In addition to phosphorylation, diffusion along the membrane is also a regulatory layer, since in-plane diffusion requires proteins to adhere to the membrane through specialized protein domains. The model, therefore, considers diffusion coefficients, which may vary depending on the membrane domain, and the dwell time of proteins on the membrane as parameters, since not all proteins are anchored in the same way, as shown in [60], for example.

At the nanoscale, as seen in Fig. 1, adhesive structures are segregated into distinct functional nanoclusters or nano-aggregates, some of which are elementary units where a protein can be phosphorylated on specific residues. From now on, we refer to a nanocluster as the elementary functional unit where a protein can be phosphorylated at some rate $\mu_p$. In a numerical analysis, where the two-dimensional space is partitioned into units, this unit is referred to as an elementary cell. Each elementary cell models a protein kinase nanocluster, and the set of elementary cells is the activation disk, loosely representing adhesive structures. The size of the activation disk, or the number of elementary units, is a parameter of the model, and their number ranges from a few to dozens. In short, the disk is a functional model for a signaling platform with a
group of immobilized kinases.

Proteins enter and exit by 2D diffusion, and they can be (de)phosphorylated at some rate ($\mu_p$ and $\mu_{dp}$), as shown in Fig. 2. These rates are parameters of the model and are proportional to the catalytic rates in the non-saturated regime [34]. They can also vary over time due to a fast-varying control variable, such as a change in conformation. This is the dynamic disorder case, a problem tackled here from a numerical perspective. For analytical calculations, we will use the effective medium approximation, which here consists of averaging the phosphorylation rate over all cells. This approximation holds if $\mu_p$ fluctuates sufficiently rapidly, and we show how this process limits the accuracy of this approximation in the probability calculation.

After the start signal at $t = 0$, the process ends in two cases: either the protein is phosphorylated, or it desorbs from the membrane at a rate $\mu_d$. What is the probability of phosphorylation? What is the average phosphorylation time? Note that if we allow loop trips by diffusion between the disk and its outside, the random walk being recurrent in 2D [15], we must impose a finite rate of desorption from the membrane to the cytosol; otherwise, this probability is 1. For a finite desorption rate, the protein is said to be mortal [61].

To better understand the role of the parameters, it is useful to first consider i) a constant diffusion constant $D$ along the membrane and ii) the two limiting cases for the completion time probability:

1. For immortal proteins (model A): We impose absorbing boundary conditions with zero probability at the disk boundary $r = R_a$. The process is complete and stops when the protein is phosphorylated (success) or when it reaches the boundary of the disk (failure). The problem is characterized by two length scales: the activation length, $l_p = \sqrt{D/\mu_p}$, and $R_a$. Since the solution to the problem involves dimensionless factors, the ratio $l_p/R_a$ is the only parameter. To within a numerical factor, this ratio is the square root of the diffusion time over the average phosphorylation time at a single site.

2. For mortal proteins (model B): The process comes to an end when either the protein is phosphorylated (success) or it desorbs from the membrane (failure). The desorption length $l_d = \sqrt{D/\mu_d}$ is the third length scale. We assume $l_d > R_a$. A mortal protein can perform many round trips before being phosphorylated if its desorption rate is sufficiently small. This is a loop or round-trip effect.

By analogy with a biological context, the differences between the two models can be visualized as follows. The first corresponds to a protein recruited from the cytosol that can bind other elements in the signaling platform (as observed for proteins able to both bind other proteins and the 2D slowly diffusing membrane of the platform). The second corresponds to the more general situation where the protein can return many times to the activation disk before desorbing into the cytosol.

Until now, we have assumed equal, constant diffusion coefficients inside and outside the activation disk. This approximation is questionable since some proteins can interact with

multiple partners before finding their target following phosphorylation. The problem is generally solved using an effective diffusion constant. This point is considered in a third section, where we show that using different diffusion coefficients $D_{1,2}$ for the inside and outside of the disk enables a crossover between the two limiting cases of models A and B above.

<figure><img src="image_2.png" /><figcaption>FIG. 1. An activation disk is a region divided into elementary cells where a protein can diffuse with diffusion coefficient $D = D_1$ and be phosphorylated at a rate of $\mu_p$. Each cell can oscillate over time between two states: active ($\mu_p > 0$), red (light gray) symbols, and inactive ($\mu_p = 0$), green (dark gray) symbols. This description is suitable for stochastic calculations. For the analytical approach, the effective medium approach means averaging the phosphorylation rate over the entire activation disk (yellow disk). Dephosphorylation will be introduced in a similar way with a local dephosphorylation rate $\mu_{dp}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Four elementary processes for a protein initially deposited in the activation disk (yellow, see Fig. 1) with diffusion coefficient $D_1$ inside the disk and $D_2$ outside: (1) The protein is phosphorylated before exiting the disk; (2) The protein is phosphorylated after a round trip (loop effect); (3) The protein is desorbed from the membrane outside the disk. Event (4) with the outside dashed circle only concerns numerical simulations and can be set aside on first reading, see Appendix A.</figcaption></figure>
### III. THE MONO-PHOSPHORYLATION CASE

<figure><img src="image_4.png" /><figcaption>FIG. 3. Probability of a protein being phosphorylated once as a function of the parameter $l_p/R_a$, where $l_p/R_a \ll 1$ corresponds to the strong phosphorylation regime. The bottom curve corresponds to model A (immortal protein) with absorbing boundary conditions at the boundary $r = R_a$, see (4). Points correspond to numerical data with different values of $D$, $\mu_p$ and $R_a$ to illustrate how different systems collapse on the same curve. The two upper curves correspond to mortal proteins (model B, see (7)) for two values of $R_a = 3$ (upper curve, $l_d/R_a \approx 1.89$) and $R_a = 6$ (median curve, $l_d/R_a \approx 0.94$). The smaller the radius $R_a$, the greater the effect of round trips (loop effect).</figcaption></figure>

In both cases, immortal and mortal proteins, not all proteins are phosphorylated after random deposition in the disk. In the mono-phosphorylation case, the calculation amounts to determining the propagator of a diffusion equation

$$ \frac{\partial \mathcal{P}_1}{\partial t} = D \Delta \mathcal{P}_1 - \mu_{1,2} \mathcal{P}_1(r,t) \quad (1) $$

with the initial condition

$$ \mathcal{P}_1(r,t=0) = \frac{1}{2\pi r_0} \delta(r-r_0) \quad (2) $$

and where the subscript 1 refers to the inside of the disk, $\mu_1 = \mu_d + \mu_p \approx \mu_p$, and 2 to the outside, $\mu_2 = \mu_d$. To calculate the probability $\mathcal{P}_1$ to be phosphorylated once, one calculates

$$ \mathcal{P}_1 = \frac{2\mu_p}{R_a^2} \int_0^{R_a} dr_0 r_0 \int_0^{\infty} dt \mathcal{P}_1(r_0, t) \quad (3) $$

to average over the starting point $0 \le r_0 \le R_a$.

#### 1. Model A

Using the calculations in Appendix B, it is first interesting to discuss the case of an immortal protein with absorbing boundary conditions (model A). As pointed out before in other contexts [62, 63], this probability is a function of only the dimensionless ratio $l_p/R_a$

$$ \mathcal{P}_1(A) = 1 - 2 \frac{l_p I_1(R_a/l_p)}{R_a I_0(R_a/l_p)} \quad (4) $$

where $I_{0,1}$ are modified Bessel functions, as shown in Fig. 3, and where the strong phosphorylation regime corresponds to small values of the parameter $l_p/R_a$.

This scaling can be made more explicit by noting that the ratio $l_p/R_a = \sqrt{D/(\mu_p R_a^2)}$ is, within some numerical factor, the square root of the fraction of time it takes for the protein to be phosphorylated without outflow to the typical time it takes for a free protein to diffuse out of the activation disk. Let $N_a$ be the total number of kinase proteins in the disk. By definition, $\mu_p$ is an average rate and $\mu_p = N_a k_a h^2 / R_a^2$, where $k_a$ is the production rate per enzyme of the product in the Michaelis-Menten kinetics [34] and $h$ the size of a functional enzyme. Thus $l_p/R_a$ depends on $N_a$ but not on $R_a$. As a consequence, the scaling plot in Fig. 3 shows that the probability for phosphorylation is independent of the size of the activating disk for a given number of kinases, because increasing the concentration by decreasing the size does not change the ratio of these two time scales.

However, because diffusion takes time, the size of the disk does matter for the mean phosphorylation time $<t_{phos}>$, and it increases with the size $R_a$. For future reference, the natural scale for this process is the mean escape time through diffusion from the disk, $t_{diff} = R_a^2/8D$, where the factor 8 arises because one averages over the starting point.

To calculate this time, we refer to the calculation in Appendix C. The probability $\mathcal{P}_1(A,t)dt$ of the protein being phosphorylated once between $t$ and $t+dt$ is

$$ \mathcal{P}_1(A,t)dt = \frac{1}{\mathcal{Z}} \mu_p e^{-\mu_p t} (1 - \mathcal{F}_{\mu_p=0}(t))dt \quad (5) $$

In (5), $\mathcal{F}_{\mu_p=0}(t)$ is the probability to cross the boundary of the disk for the first time at exactly $t$, and $\mathcal{Z}$ is a normalization factor.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Average time for a protein to be phosphorylated at a single site, normalized by diffusion time, see (C13). As indicated by the points obtained numerically, this curve, once normalized, is a scaling curve as a function of $l_p/R_a$. The dotted curve corresponds to $<t_{phos}> = 1/\mu_p$, which would be the value of the mean phosphorylation time for a disk of infinite radius. Down triangles and circles correspond to different sizes $R_a$ and make the scaling property explicit. For $l_p/R_a > 1$, the probability that a protein will be phosphorylated is negligible.</figcaption></figure>

Taking the average value over the subset of proteins that
<figure><img src="image_6.png" /><figcaption>FIG. 5. Normalized histograms of the number of times a protein crossed the boundary of the adhesive disk before being phosphorylated. There are therefore only even values. The two plots correspond to $l_p/R_a = 0.21$ (orange - light gray) and $0.42$ (green - dark gray). For small $l_p/R_a$, in the high phosphorylation rate regime, proteins are rapidly phosphorylated. Data corresponds to $\mu_d = 5 \times 10^{-4}$, $R_a = 3$ so that $l_d/R_a = 1.89$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Dynamic disorder case: probability of a protein being monophosphorylated (model B) when the activation zone in Fig. 1 contains only 16 cells dynamically alternating between two states, one of which can phosphorylate and the other not. For this small cluster, the activation zone has the shape of a square, as seen in the insets which represent snapshots of the two configurations in the activation zone. The stochastic simulation runs with equal rates $k_+ = k_-$ between the two states. The points with symbols • correspond to slow rates $(0.01s^{-1})$, while the dashed curve corresponds to (7) with an effective rate $\mu_p/2$.</figcaption></figure>

have been phosphorylated exactly once, one finds

$$ \langle t_{phos} \rangle = t_{diff} \mathcal{T}(l_p/R_a) \quad (6) $$

where the function $\mathcal{T}(l_p/R_a)$ is given in Appendix C, see (C13), and is plotted in Fig. 4. This calculation shows that $\langle t_{phos} \rangle$ is generally smaller than $1/\mu_p$, with equality only in the very high phosphorylation rate regime. This trend can be explained by the fact that diffusion limits phosphorylation over short periods of time and skews the sampling to calculate the average by imposing a cut-off for the long times.

## 2. Model B: Homogeneous case

The specific case of model A is in stark contrast with the general result where we allow loop effects (model B) for mortal proteins ($D_1 = D_2 = D$). Referring again to Appendix B, the general result reads as

$$ \mathcal{P}_1(B) = 1 - 2 \frac{l_1 I_1(R_a/l_1) K_1(R_a/l_d)}{l_d R_a Z_1} \quad (7) $$

where $Z_1$ is a combination of Bessel functions and $l_1 = \sqrt{D/(\mu_d + \mu_p)} \approx l_p$, $\mu_d \ll \mu_p$. As it should, result (4) follows from (7) in the limit of a high desorption rate outside the activating disk, $l_d \to 0$. As shown in Fig. 3, this probability increases markedly, even for moderate values of desorption length, indicating the importance of round-trip effects. This increase in probability naturally depends on the average time that must elapse before the first phosphorylation occurs. This time, in turn, depends on the desorption rate. In the case of model B, the reference time for this process is set by the average time $t_{dwell}$ that the protein has spent in the disk when diffusing back and forth. Calculations in Appendix D show that $t_{dwell}$ is given by

$$ t_{dwell}/t_{diff} = \frac{l_d^2}{8R_a^2} \left[ 1 - 2 \frac{R_a}{l_d} I_1(R_a/l_d) K_1(R_a/l_d) \right] \quad (8) $$

with rapid variations when increasing $l_p/R_a$. A simple way to quantify this loop effect is to count the number of times a protein needs to cross the boundaries of the disk back and forth before being phosphorylated. The result is shown in Fig. 5 in the high and moderate phosphorylation regime, with the conclusion that even for moderate ratios $l_d/R_a$, the effect of allowing loops is strong due to the long tails common in these kinds of problems. The importance of these tails is particularly apparent in Fig. 5, which shows how a protein can cross the disk boundary some ten times, even though the desorption length is only twice the size of the disk. This increase in probability is naturally at the expense of the average phosphorylation time, which is greater than in the previous case.

## 3. Model B: Dynamic disorder case

Fig. 5 is obtained using a stochastic solver for the Master equation equivalent to the problem [64] and Appendix A. The advantage of the numerical approach is the ability to compare the previous results with the dynamic disorder case, where the disk has a very limited number of kinases that alternate with rates $k_+$, $k_-$ between an active state, capable of phosphorylating the protein with $\mu_p > 0$, and an inactive state, unable to do so with $\mu_p = 0$. In other words, we assume that the enzymatic activity per cell depends on a time-varying control parameter $C$, such as a change in molecular conformation, as in the classical gated binding problem [65–67]. This telegraphic noise, often called dichotomous noise, breaks detailed balance and is often used to describe non-equilibrium states[68, 69].
<figure><img src="image_8.png" /><figcaption>FIG. 7. Probability that a protein will be mono-phosphorylated in model B as a function of the average time required for the process to reach completion (success). The top x-axis corresponds to normalization by the diffusion time $t_{diff} = R_a^2/8D$, while the bottom scale is a function of the time actually spent in the nanocluster before desorption, see Eq. (D3). Due to the non-zero desorption rate, $\mu_d > 0$, not all proteins are phosphorylated. The characteristic length $l_d = \sqrt{D/\mu_d}$ is therefore finite ($D_1 = D_2 = D$). The solid line corresponds to the analytical formula (E2) and (E4) calculated to the first order in $1/\mu_p$ when $\mu_p$ is large, and it is a priori valid only for $<t_{phos}>/t_{diff} \ll 1$ (see Appendix E). The symbols • correspond to the results of a numerical simulation where the activation zone is circular with size $R_a$ and $R_a/l_d = 0.56$. For the other symbols, the activation zone has the shape of a small square of 16 elementary cells. The symbols ♦ correspond to a homogeneous system with a constant phosphorylation rate $\mu_p$, while the symbols ▲ correspond to the case where cells dynamically alternate between an active state and an inactive state (dynamic disorder).</figcaption></figure>

For fast time-varying control parameters, we expect renormalized effective rate constants. Fig. 6 plots the probability $P_1(B)$, where $l_p/R_a$ is now a function of the time-averaged enzymatic activity defined from the fraction of time spent in the active conformation:

$$ \mu_p \rightarrow \langle \mu_p \rangle = k_+ / (k_+ + k_-) \mu_p \quad (9) $$

If the fluctuations are sufficiently fast, it is sufficient to take the time average $\langle \mu_p \rangle$, and the probability coincides with (7). In the other limit of static disorder for slow variations between the active and inactive state, the two curves differ slightly, because for a low number of kinases, there is a non-zero probability that the protein will desorb without being phosphorylated, even within the limits of very high phosphorylation rates.

As mentioned earlier, an ideal single-molecule experiment measures the probability of mono-phosphorylation, $P_1(B)$, as well as the average time $\langle t_{mean} \rangle$ its take for a protein to be phosphorylated. Each of these two quantities depends on the average phosphorylation rate $\langle \mu_p \rangle$, which is one of the parameters of the problem. Consequently, the parametric curve

$$ \langle t_{mean} \rangle, P_1(B) \quad (10) $$

makes it possible to avoid determining this parameter in making predictions. As seen in Fig. 7, this graph gives an almost

straight line we now interpret.

As a reference, Appendix E provides the leading behavior of both quantities in the high phosphorylation regime. Each term varies as $1/\mu_p$, so that the graph shows a straight line consistent with the stochastic numerical data for sufficiently large radii $R_a$, even in the low phosphorylation regime. This graph shows, however, small deviations for small, dynamically disordered activation disks with a square shape (see insets of Fig. 6) only in the low phosphorylation regime. As for the problem in [62], one expects, therefore, that the asymptotic forms are independent of the circular geometry considered here. As a consequence of this parametric plot, changing the effective phosphorylation rate, for example, by modifying the residence time of kinases, allows the graph to be traversed, assuming that the desorption rate $\mu_d$ of the protein remains unchanged.

A second consequence of the approach proposed here is to show how phosphorylation clusters enable or prevent crosstalk between different signaling pathways by renormalizing effective rate constants. In parallel with single-molecule kinetics [70, 71], the inverse of the first moment of the waiting time distribution for a single protein gives an effective cluster-size-dependent phosphorylation rate for an ensemble average measurement. In the large phosphorylation regime with $\mu_p \gg \mu_d$, this rate, $\mu_{p,eff}$, is proportional to $\mu_p$ renormalized by a function of $R_a/l_d$ with

$$ \mu_{p,eff} \approx 2\mu_p \frac{K_0^2(R_a/l_d)}{Q(R_a/l_d)} \quad (11) $$

as given in Appendix E. This rescaling takes the value 1 when $R/l_d$ is large but approaches 0 in the other limit and has a simple physical interpretation. When the desorption rate is large, $R/l_d \gg 1$, phosphorylation takes place on a single site without diffusion. The rate is therefore $\mu_p$. The other limit, equivalent to that of a large diffusion coefficient, implies that the protein, although returning many times to the activation disk before desorbing, spends very little time in this disk and, as a result, its phosphorylation rate decreases.

4. Crossover between model A and B

So far, we have considered the special case of a homogeneous diffusion coefficient. In the context of cell adhesion, this assumption is not necessarily always valid.

The above results can be generalized to the following situation: For $0 \le r \le R_a$, let $D_1$ be the diffusion coefficient, while it is $D_2$ outside, with $D_2 > D_1$. Clearly, excursions away from the activation disk are favored, and the mortal protein is expected to desorb more frequently before returning in its round trips. Using the results of Appendix B 2, Fig. 8 clearly shows the effect of increasing $D_2$ relative to $D_1$, where the probability of phosphorylation crosses over between the two limiting cases, model A and model B.
<figure><img src="image_9.png" /><figcaption>FIG. 8. Crossover between model A and B. Probability of a protein being phosphorylated once as a function of the parameter $l_p/R_a$, with $l_p = \sqrt{D_1/\mu_p}$. For reference, the upper dashed curve corresponds to $D_2 = D_1$ (model B) and the continuous bottom one is for absorbing boundary conditions (model A) at the boundary of the activation disk. The two other dashed curves correspond to $D_2 = 10D_1$ and $D_2 = 100D_1$ (bottom) and illustrate how the system crosses over from model A to B.</figcaption></figure>

## IV. THE MULTI-PHOSPHORYLATION CASE

We consider a protein that can be phosphorylated on an arbitrary number of residues and take $D_2 = D_1 = D$. Dephosphorylation must be taken into account and will be addressed in a second subsection. We introduce a threshold $n$ as the minimum number of phosphorylated residues to serve as a criterion for process completion. If a protein is deposited at time $t = 0$ in the disk, the process is completed in two cases:

1. Model A (immortal protein): The protein exits the disk without having been phosphorylated on $n$ residues (failure), or the protein is phosphorylated on $n$ residues without having crossed the boundary (success).

2. Model B (mortal protein): The protein is phosphorylated on $n$ residues (success) or desorbs from the membrane (failure).

### A. Simple multi-phosphorylation networks

Whatever the model (A or B), multi-phosphorylation can be either sequential or random. In the former case, the sequence of events follows an order where, for example, residue 1 must be phosphorylated for residue 2 to be phosphorylated, and so on. In an alternative version, the order in which the residues are phosphorylated is irrelevant, and the networks are said to be random, as shown in Fig. 9.

Each network behaves differently, since the phosphorylation propensity is constant in the sequential case, whereas it varies at each level of the hierarchy for the random network. For example, in the random case, to go from level $k = 0$ (no phospho-residues) to level $k = 1$ (1 phosphorylated residue among $n$), this propensity is $n\mu_p$. This difference

translates into a higher probability for the random network with a shorter time, given that the process is completed as soon as $n$ phospho-residues are reached (success). This is illustrated in Fig. 10 for a mortal protein that can exit and re-enter the activation disk before desorption (failure). Thus, a sequential network has a lower probability of being fully phosphorylated, and this takes a longer time than for a random process.

To calculate this probability, let us first define the propensity function [58] to go from level $n-k$ to level $n-k+1$ of the hierarchy. The probability of a site being phosphorylated per unit time is $\mu_p$. So the probability of one of the residues being phosphorylated per unit time is the sum of these probabilities, i.e., $(n-k)\mu_p$, and the total propensity is

$$P(n-k \rightarrow n-k+1) = \mathcal{P}_1(k\mu_p) \neq k\mathcal{P}_1(\mu_p) \quad (12)$$

Consequently, multi-phosphorylation is not additive, and because $\mathcal{P}_1$ depends on the diffusion constant $D$, diffusion induces some cooperative effects between the residues. For simplicity, consider model A. Let $\mathcal{P}_1(\mu_p)$ be the probability (before escape) if it has only one possible phosphorylation site. Then, as all residues play the same role, the protein has a probability $n \times (1/n) \mathcal{P}_1(n\mu_p)$ of being phosphorylated at at least 1 site before exiting. The probability to exit the disk with exactly one site among the $n$ is therefore

$$\mathcal{P}_1(n\mu_p)(1 - \mathcal{P}_1((n-1)\mu_p)) \quad (13)$$

with a probability $1 - \mathcal{P}_1(n\mu_p)$ of escaping without being phosphorylated. The process is the same at each step of the hierarchy $k = 1, 2, \dots, n$ phosphorylated residues. The probability of escaping with exactly two phosphorylated residues is

$$\mathcal{P}_1(n\mu_p)\mathcal{P}_1((n-1)\mu_p)(1 - \mathcal{P}_1((n-2)\mu_p)) \quad (14)$$

and for exactly all $n$ phosphorylated residues

$$\mathcal{P}_n = \mathcal{P}_1(n\mu_p)\mathcal{P}_1((n-1)\mu_p)\dots\mathcal{P}_1((n-n+1)\mu_p) \quad (15)$$

We verify

$$\begin{aligned} & 1 - \mathcal{P}_1(n\mu_p) + \mathcal{P}_1(n\mu_p)(1 - \mathcal{P}_1((n-1)\mu_p)) \\ & + \mathcal{P}_1(n\mu_p)\mathcal{P}_1((n-1)\mu_p)(1 - \mathcal{P}_1((n-2)\mu_p)) \\ & + \dots + \mathcal{P}_1(n\mu_p)\mathcal{P}_1((n-1)\mu_p)\dots\mathcal{P}_1(1 \times \mu_p) \\ & = 1 \end{aligned} \quad (16)$$

Using previous results for $\mathcal{P}_1$, we can calculate the probability that the protein is phosphorylated on $k = 0, 1, 2, \dots, n$ residues. This is shown in Fig. 11 for model A and in Fig. 10 (inset) for model B. In both cases, the system behaves as an all-or-nothing switch as the scaling parameter $l_p/R_a$ is varied. At a sufficiently high rate of phosphorylation, either all of the $n = 10$ residues are phosphorylated, or none of them are. The size $R_a$ of the nanoclusters cannot select the number of phosphorylated residues. Because the propensity for phosphorylation varies in the hierarchy, the time necessary for the process to go from level $k = 0$ to full phosphorylation with $k = n$ is
<figure><img src="image_10.png" /><figcaption>FIG. 9. Phosphorylation networks. A protein can be phosphorylated on any or all of its *n* residues (*n* = 3 here) as indicated by the arrow. The left-hand column shows the total number *k* of phosphorylated residues at each level of the hierarchy. The network in the middle is a sequential network, while the network on the right is a random network. For the random network, there are $C_n^k$ possibilities at each level for a total of 2ⁿ possible states and 2ⁿ⁻¹*n* phosphorylation reactions. Dephosphorylation can be added either sequentially by modifying each edge (arrow) to make it bidirectional, e.g., (001) ↔ (011), or by adding a complete dephosphorylation for each residue, for example, (011) → (000), both models with a $\mu_{dp}$ rate.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Mean time for a protein to be phosphorylated on *n* = 10 residues (normalized by $t_{\mu_p=0}^{in}$ as defined in (D3)). The upper curve (circles) corresponds to the sequential network and the bottom one to the random network (up-triangles). Inset: Probability for a protein to be phosphorylated on *n* = 10 residues for a sequential (circle, bottom curve) and random (up-triangle, upper curve) network. Numerical data fit the two analytical curves calculated in the text with error due to finite-size effects (see (15)). Data correspond to mortal proteins with $\mu_d = 5 \times 10^{-4}$ or $l_d/R_a = 1.89$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11. Probability of exiting the disk for the first time with $k=1,...,5$ or 10 phospho-residues ($R=3, \mu_d=0.0$, immortal proteins). The phosphorylation network is random as in Fig. 6 ($n=10$). For high phosphorylation rates, $l_p/R \ll 1$, all residues are phosphorylated before crossing the boundary (see $k=10$). Inset: same probabilities from stochastic simulations ($k=10$: red (dark gray) circles, $k=0$: +, $k=1$: ×, etc.).</figcaption></figure>

much faster in the random multi-phosphorylation case than in the sequential model, as shown in Fig. 10.

As before and to compare with Fig. 7, Fig. 12 shows the graph of the probability that a protein is phosphorylated on *n* = 10 residues as a function of the average time to complete the process. Clearly, the linear nature of the behavior is still present within the limits of the phosphorylation rates. However, as expected, the characteristic times for the kinetics are significantly slower in the case of multi-phosphorylation.

## B. Multi-phosphorylation with phosphatases

Dephosphorylation, a function orchestrated by phosphatases, can be easily included in the model. The first way to do this is to assume that phosphorylation sites in the disk $R_a$ are also dephosphorylation sites, so that a protein visiting the $R_a$ disk by diffusion is subject to a chemical reaction of the type

$$ (\dots, 0_i, \dots) \xrightarrow{\mu_p} (\dots, 1_i, \dots) \xrightarrow{\mu_{dp}} (\dots, 0_i, \dots) \quad (17) $$

with rates $\mu_p$ and $\mu_{dp}$. Here, this scheme is termed sequential, because each phospho-reaction has its own inverse. To parallel the McKeithan proofreading mechanism, another scheme, referred to here as total, includes complete dephosphorylation
<figure><img src="image_13.png" /><figcaption>FIG. 12. Probability that a protein is phosphorylated on n = 10 residues as a function of the mean time necessary to reach completion (success). As in Fig. 7 and for comparison, the time scale is normalized in two ways. The top x-axis corresponds to normalization by the diffusion time $t_{diff} = R_a^2/8D$, while the bottom scale is a function of the time actually spent in the nanocluster before desorption, see Eq. (D3). Data correspond to $R_a/l_d = 1.78$ with $R_a = 3$.</figcaption></figure>

at all levels of the hierarchy:

$$ (\dots, 1_i, \dots, 1_j, \dots) \xrightarrow{\mu_{dp}} (0_1, \dots, 0_i, \dots, 0_j, \dots, 0_n) \quad (18) $$

It is instructive to start with the first scheme by generalizing model A. The process is completed if the number of phosphorylated residues reaches a maximum threshold, or if the protein crosses the boundary. In the latter case, the number of phospho-residues is retained. Using (17), we can plot the histograms shown in Fig. 13, where all results correspond to the same phosphorylation rate $\mu_p$. For this value of $\mu_p$, more than 90% of proteins are fully phosphorylated if dephosphorylation is not included in the model. Clearly, dephosphorylation, the propensity of which increases with the number of phosphorylated residues, enables histograms to be shifted leftward, and therefore the disk determines a specific number of residues to be phosphorylated in the hierarchy $k = 0, \dots, n$. In this case, the size $R_a$ is a selection factor, but as seen below, this case is marginal.

Taking into account the loop effects of model B changes this picture. For the same dephosphorylation scheme (17), Fig. 14 shows the effect of sequential dephosphorylation for a protein to be completely phosphorylated. In this case, dephosphorylation mainly affects the subset that goes back and forth, since for $l_p/R_a \ll 1$, proteins are phosphorylated very rapidly and the tangent to the curve remains horizontal. Increasing $\mu_{dp}$ very quickly, the curves take on the shape of an inverted sigmoid, with a threshold effect for a value of parameter $l_p/R_a$ below which proteins are completely phosphorylated. The effect is the same for the total phosphorylation scheme of Eq. (18), see Fig. 15, and is almost independent of the maximum number of residues, since decreasing the maximum number of residues from $n = 10$ to $n = 5$ produces the same effect, as seen in Fig. 16. We therefore conclude that this is a robust property for networks.

To highlight the importance of the local nature of post-

<figure><img src="image_14.png" /><figcaption>FIG. 13. Normalized histograms of the number n of phosphorylated residues when the protein first crosses the disk boundary ($n = 10$). The sequential dephosphorylation rate is set at $\mu_{dp} = 0.05$ with three phosphorylation rates ($\mu_p = 0.25$ (green - dark gray), $\mu_p = 0.05$ (orange - medium gray) and $\mu_p = 0.01$ (yellow - light gray)), corresponding to $l_p/R_a = 0.084, 0.18, 0.42$. In this model, the dephosphorylation rate selects the number of phosphorylated residues, since the histograms are either centered on small n's ($\mu_p \ll \mu_{dp}$) or on large n's ($\mu_p \gg \mu_{dp}$).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 14. Probability for $n = 10$ residues to be phosphorylated for different sequential dephosphorylation rates ($\mu_{dp} = 0.00, 0.01, 0.025, 0.05, 0.1$ and $\mu_d = 5 \times 10^{-4}$). The line corresponds to no dephosphorylation rate $\mu_{dp} = 0$ and serves as a reference. The higher the rate of dephosphorylation, the more the curve is sigmoid-shaped. Inset: Probability of the residue number k being phosphorylated ($\mu_{dp} = 0$). In the high phosphorylation limit, $l_p/R_a = 1$, the protein is fully phosphorylated on $k = 10$ residues (blue down triangles). In the other limit, the protein is moderately phosphorylated at all of the residues ($k = 0$, symbol +).</figcaption></figure>

translational modifications in obtaining switch-like responses, consider Fig. 17. In this case, we depart from the previous rule and introduce phosphatases only outside the activation circle, so that only proteins localized outside the disk can be dephosphorylated in a single reaction (total dephosphorylation). Even for very high rates of dephosphorylation, the response is gradual, not abrupt, because dephosphorylation only affects the tail end of the distribution where proteins are most subject to round-tripping.
<figure><img src="image_16.png" /></figure>

FIG. 15. Probability for $n = 10$ residues to be phosphorylated for different total dephosphorylation rates ($\mu_{tdp} = 0.0, 0.01, 0.05, 0.1$ (magenta circles - dark gray)), see Eq. (18). The line corresponds to no dephosphorylation rate $\mu_{tdp} = 0$ with yellow (light gray) lozenges and serves as a reference. As in Fig. 14 for sequential dephosphorylation, the higher the rate of dephosphorylation, the greater the variations in the curve, so that the system switches abruptly from full phosphorylation to zero by varying the parameter $l_p/R_a$.

<figure><img src="image_17.png" /></figure>

FIG. 17. Probability of phosphorylating residue k for a protein with a maximum number of n = 10 residues (k = 10 corresponds to ▼ and k = 0 to + symbols). For this simulation, the protein can only be (multi)phosphorylated in the activation disk of radius $R_a$, whereas it can only be dephosphorylated (total dephosphorylation) outside it (area colored red). The rate of $\mu_{tdp}$ dephosphorylation is 5 times greater than the maximum $\mu_{dp}$ in Fig. 14 and only the tail of the distribution without dephosphorylation (solid curve) is affected by dephosphorylation.

<figure><img src="image_18.png" /></figure>

FIG. 16. Probability for $n = 5$ residues all to be phosphorylated (upper green losange symbols, the maximum number of residues is $n = 5$). For reference, this probability, but for a maximum number of $n = 10$ residues also all phosphorylated, is shown as blue circle symbols (same as in the inset of Fig. 10). The orange pentagon and magenta triangle symbols correspond to a dephosphorylation rate $\mu_{dp} = 0.1, 1$ for a maximum number of $n = 5$ also all phosphorylated.

V. POTENTIAL BIOLOGICAL RELEVANCE AND CONCLUSION

The ability to mono- or multi-phosphorylate proteins is one of the essential features of cell signaling networks, and this work proposes a reaction-diffusion approach to study this problem in the context of adhesion. To calculate these probabilities, we have introduced a dimensionless parameter that relates individual properties of each protein, such as the diffusion coefficient, to the spatiotemporal average enzymatic activity in an activation zone, as well as to its characteristic size. This probability is highly dependent on this parame-

ter and, depending on phosphatase activity, can even lead to abrupt or almost switch-like biochemical behavior, in contrast to more gradual variations. Introducing this effective parame- ter together with characteristic length and time scales is, in our view, a way of rationalizing the processes that control molecular self-assembly. It allows also to define effective nanoscopic rates to rationalize crosstalks between different signaling pathways.

It is likely that multi-phosphorylation plays a different role than mono-phosphorylation in the stability of molecular as- semblies and biomolecular condensates that the cell regulates, among others, in the immediate vicinity of the adhesive zones. In analogy to what was shown for the SH3-Proline-rich region [72], phase separation following multi-phosphorylation could be driven by the assembly of multivalent proteins. Multi- phosphorylation of residues, larger than 3-4, will support a tunable multivalency on a protein of interest that is essential for phase stability. Moreover, this phosphorylation-dependent multivalency could also provide greater local molecular rigid- ity and thus be implicated in a new type of mechano-sensing. In terms of mechanobiology, control of multi-phosphorylation is essential to activate a specific signaling pathway through p130Cas [73].

This work demonstrates that simple principle models can account for a fine-tuned regulation of mono- versus multi- phosphorylation, taking into account the effect of the envi- ronment and involving biochemical switches working with a small number of molecules. From a biochemical point of view, the main hypothesis is that phosphorylation is a locally controlled phenomenon, as for Tyrosine residues, whereas it seems to be more diffuse on a cellular scale for Ser-Thr ones because of specificity or abundance reasons [74].

Controlling the enzymatic activity within adhesive struc- tures is obviously the nerve center of any biochemical control,
and $\mu_p$ is interpreted here as an effective rate constant that may depend on a "control variable," say $\mathcal{C}$, such as a change of conformation or a local concentration of co-activators. Triggering a signaling pathway implies, however, some kind of time averaging. As usual [65–67], if $\mathcal{C}(t)$ varies in time and if its fluctuations are sufficiently fast with respect to the diffusion time, we have shown using spatiotemporal variations of $\mu_p(\mathcal{C})$ among elementary cells that $\mu_p$ can be replaced by the time-averaged effective rate, which reflects the dynamics of enzymatic activity. Changes in this dynamics by adding telegraphic (dichotomous) noise can have many origins, for example, changes in the fraction of time an enzyme spends in its active conformation or changes in how proteins are recruited. Since most of the kinases are cytosolic proteins recruited from the cytosol, Src-kinase being a paradigmatic example, changing their mode of recruitment by molecular engineering will bias their dynamics with a potential strong effect on their effective enzymatic activity. Thus, our model, which demonstrates abrupt variations in phosphorylation, represents a first step towards a potential explanation of the major changes in molecular assemblies demonstrated in [21] where opto-Src lacks its major anchor domain [75, 76]. Note that this approach, which focuses on the dynamics and the labile character of a single enzyme to induce changes in its function, has also been proposed in a parallel context in [23] for a GTP-ase.

Triggering a signaling pathway also implies some kind of local space averaging, and the second principle on which we have based our analysis is that this sensing specificity occurs through effective parameters to allow crosstalk. One signaling pathway interferes with another, as it reads biophysical parameters influenced by the other. The integration of their relative dynamics parameters should control their level of crosstalk and thus drive the system from model A to B with a switch-like behavior. The third principle highlighted in this work is time. To test proteins with nanoclusters having a short lifetime, multi-phosphorylation must occur rapidly, with characteristic times that depend on the ensemble of phosphorylable substrates available.

Following on from works on signal transduction [3, 77], this work suggests that it is possible to induce two different signaling events by playing on the spatial distribution and spacing of phosphorylation sites on a single substrate, leading to structural and temporal changes in the downstream molecular complexes associated with these substrates.

As studied here, however, cooperative effects that translate into parameter-dependent abrupt changes, enabling fine regulation by multi-phosphorylation, must include phosphatase networks, whose effects are less often discussed in *in vivo* studies than those of kinases.

## ACKNOWLEDGMENTS

O.D. is partly supported by the Agence Nationale de la Recherche PRC program (ANR-22-CE13-0046). We thank G. Gogl for useful discussions and an anonymous referee to draw our attention to the problem $D_2 \neq D_1$. We would also like to

thank the second referee for his comments, which helped us improve the content of this paper, and for drawing our attention to Ref. [30].

## Appendix A: Stochastic simulations

The stochastic counterpart to the Fokker-Planck equation is the Master equation, which we solve numerically using the Gillespie SSA algorithm. We follow the approach of [78] with an updated list of references from [79, 80]. The method is also described in [64] and is a spatial stochastic simulator for chemical reaction networks, see also Ref. [81].

The computational domain is divided into a $25 \times 25$ grid of compartments to simulate diffusion as Markov jumps between compartments and reactions within each compartment. Each compartment is treated as spatially homogeneous, but concentrations can differ between them. If $A_i$ is the number of molecules of chemical species A in the i-th compartment, the diffusion-reaction process is described by the following system of chemical reactions:

$$A_1 \underset{D/h^2}{\overset{D/h^2}{\rightleftharpoons}} A_2 \underset{D/h^2}{\overset{D/h^2}{\rightleftharpoons}} A_3 \dots A_{N-1} \underset{D/h^2}{\overset{D/h^2}{\rightleftharpoons}} A_N \quad (A1)$$

Here, $h$ is the size of one compartment and $D$ is the diffusion constant in the continuum case. Typical simulations were run with $D/h^2 = 0.1$, which provides a reference rate for a numerical experiment [78].

At each time step, the algorithm computes a propensity matrix to determine the time and type of the next chemical reaction. To save computational time, the propensity matrix for diffusion is updated at each step using a complex double-linked list of structures. The program can process any binary reaction with local reaction constants and generates the in silico statistics used in the text. This approach handles the combination of different states directly without approximation. For instance, the case $n = 10$ requires 1025 different protein phosphorylation classes with 6144 reactions. Because this method explicitly includes protein diffusion, it naturally stands out from other well-mixed approaches. Trajectory animations can be generated at each time step using the Dislin library [82].

Since the diffusion domain in a numerical simulation is necessarily finite, we introduce an absorbing boundary on a circle with a radius larger than that of the activation disk. A counter is used to track events where the protein makes excursions wider than the numerical diffusion domain allows (event (4) in Fig. 2 on the outer circle). These events are discarded when averaging over different realizations, which introduces a controllable numerical error. The other parameters are the same as in the text, with averages taken over $20 \times 10^3$ to $10^6$ samples, depending on the case.
## Appendix B: Propagator

Here we demonstrate (4) and (7). We place a protein at $r = r_0$, $0 \le r_0 \le R_a$, at $t = 0$. Calculating the probability of the protein being phosphorylated is equivalent to calculating

$$ 
\begin{aligned}
& 2\pi\mu_p \int_0^R dr r \int_0^\infty dt \mathcal{P}_1(r,t) \\
& = 2\pi\mu_p \int_0^R dr r \tilde{\mathcal{P}}_1(r,s=0)
\end{aligned} 
\quad (B1) $$

where $\mathcal{P}_1(r,t)$ is the survival probability and $\tilde{\mathcal{P}}_1(r,s=0)$ its Laplace transform. We average the results over $r_0$ and we separate into regions 1, 2 with the definitions ($D_{i=1,2}$ are diffusion coefficients)

$$ 
\begin{aligned}
\text{Region 1: } D_1, 0 \le r \le R_a \\
\text{Region 2: } D_2, R_a < r \le \infty
\end{aligned} 
\quad (B2) $$

where the Laplace transform is solution of:

$$ \frac{1}{r} \frac{\partial}{\partial r} r \frac{\partial \tilde{\mathcal{P}}_1}{\partial r} - q_1(s)^2 \tilde{\mathcal{P}}_1 = - \frac{1}{2\pi r_0 D_1} \delta(r - r_0) \quad (B3) $$

with

$$ q_{1,2}^2(s) = \frac{\mu_{1,2} + s}{D_{1,2}} \quad (B4) $$

In the notations of the text $\mu_1 = \mu_p + \mu_d \approx \mu_p$ (phosphorylation rate, desorption rate) and $\mu_2 = \mu_d$. We have

$$ q_1(s=0) = 1/l_p \quad q_2(s=0) = 1/\lambda \quad (B5) $$

We look for a solution in the form

$$ 
\begin{aligned}
A_1 I_0(q_1(s)r) \text{ for } r < r_0 \\
A_2 I_0(q_1(s)r) + B_2 K_0(q_1(s)r) \text{ for } r_0 \le r < R_a \\
B_3 K_0(q_2(s)r) \text{ for } r \ge R_a
\end{aligned} 
\quad (B6) $$

with boundary conditions:

$$ \left. \frac{d \tilde{\mathcal{P}}_1}{dr} \right|_{r=r_0^+} - \left. \frac{d \tilde{\mathcal{P}}_1}{dr} \right|_{r=r_0^-} = - \frac{1}{2\pi D_1 r_0} \quad (B7) $$

where probabilities are continuous in $r = R_a$.

It will be useful to define:

$$ Z_1 = \frac{D_2}{D_1} q_2(s) I_0(q_1(s)R_a) K_1(q_2(s)R_a) + q_1(s) K_0(q_2(s)R_a) I_1(q_1(s)R_a) \quad (B8) $$

$$ Z_2 = q_1(s) K_0(q_2(s)R_a) K_1(q_1(s)R_a) - q_2(s) \frac{D_2}{D_1} K_0(q_1(s)R_a) K_1(q_2(s)R_a) \quad (B9) $$

Hence:

$$ 
\begin{aligned}
A_1 &= \frac{1}{2\pi D_1} K_0(q_1(s)r_0) + \frac{Z_2}{2\pi D_1 Z_1} I_0(q_1(s)r_0) \\
A_2 &= \frac{Z_2}{2\pi D_1 Z_1} I_0(q_1(s)r_0) \\
B_2 &= \frac{1}{2\pi D_1} I_0(q_1(s)r_0) \\
B_3 &= \frac{1}{2\pi D_1 R_a Z_1} I_0(q_1(s)r_0)
\end{aligned} 
\quad (B10) $$

It's easier to do the mean value on $r_0$ first:

$$ I(r) = \frac{2}{R_a^2} \int_0^{R_a} dr_0 r_0 \tilde{\mathcal{P}}_1(s=0, r) \quad (B11) $$

and then do the integral over $r$ times $\mu_p$

$$ 2\pi\mu_p \int_0^{R_a} dr r I(r) \quad (B12) $$

to find the Laplace transform at $s = 0$ of a probability of a protein being phosphorylated

$$ \mathcal{P}_1(B) = \frac{\mu_p}{q_1^2 D_1} \left[ 1 - \frac{2q_2}{q_1 R_a Z_1} \frac{D_2}{D_1} I_1(q_1 R_a) K_1(q_2 R_a) \right] \quad (B13) $$

with

$$ q_{1,2} = q_{1,2}(s=0) = \sqrt{\mu_{1,2}/D_{1,2}} \quad (B14) $$

### 1. Limit cases with equal diffusion constant $D_1 = D_2 = D$

Take $D_{1,2} = D$ first. We find the limit of problem A (absorption on the disk boundary $r = R$), when $\mu_2 = \mu_d \gg 1$, i.e. $q_2 \gg 1$. In this limit

$$ Z_1 \approx q_2 I_0(q_1 R_a) K_1(q_2 R_a), \quad q_2 \gg 1 \quad (B15) $$

which gives the desired result, since $q_1^2 = \mu_p/D = 1/l_p$.

$$ 1 - \frac{2l_p I_1(R_a/l_p)}{R_a I_0(R_a/l_p)} \quad (B16) $$

In this limit, the probability flux across the boundary of the disk gives the probability to exit for the first time at time $t$. Its Laplace transform is

$$ \tilde{\mathcal{F}}_{\mu_p}(s) = \frac{2l(s) I_1(R_a/l(s))}{R_a I_0(R_a/l(s))} \quad (B17) $$

with

$$ l(s) = \sqrt{D/(\mu_p + s)} \quad (B18) $$

In the other limiting case where $\mu_d \ll 1$ outside the disk of size $R$, $Z_1 \approx -q_1 I_1(q_1 R_a) \ln q_2 \gg 1$ and we find that the probability is 1 in the limit where the logarithm is sufficiently large.

Surprisingly, the inverse Laplace transform $\tilde{\mathcal{F}}_{\mu_p}(t)$ of $\tilde{\mathcal{F}}_{\mu_p}(s)$ does not seem to be known. For $t \to 0$, i.e. $s \to \infty$

$$ \tilde{\mathcal{F}}_{\mu_p}(t) \approx \left( \frac{4D}{\pi R_a^2} \right)^{1/2} \frac{e^{-\mu_p t}}{t^{1/2}} \quad (B19) $$

We can also calculate the average time $t_{diff}$ it takes for the protein to leave the disk for the first time when $\mu_p = 0$ by calculating $-\tilde{\mathcal{F}}'_{\mu_p}(s=0) = R_a^2/(8D)$, which differs by a factor of 2 from the definition of the non-averaged diffusion time $R_a^2 = 4Dt_d$.
## 2. Case where $D_2 > D_1$

We study the case where the diffusion constant $D_2$ outside the disk $0 \le r \le R_a$ is much larger than the one inside the disk $D_1$. In (B13), take $D_2 \gg 1$, then $q_2 \ll 1$ and $K_1(q_2R_a) \approx 1/q_2R_a$. With logarithmic corrections $Z_1 = D_2/(D_1R_a)I_0(q_1R_a)$, and $K_1(q_2R_a) \approx 1/q_2R_a$. So we find again the result of Eq. (B16). The crossover between $D_2 = D_1$ (model B) and $D_2 \gg D_1$ (model A) is illustrated in Fig. 8.

# Appendix C: Mean time (model A)

This Appendix takes up the previous reasoning on calculating the probability of phosphorylation. We deposit the protein inside the yellow disk at a distance $a < R$ from the center and we calculate the probability that this protein reaches the boundary $r = R$ between $t$ and $t+dt$. This probability is denoted $\mathcal{F}_{\mu_p}(t, r_0)dt$ and its Laplace transform is $\tilde{\mathcal{F}}_{\mu_p}(s, r_0)$. We have

$$ \tilde{\mathcal{F}}_{\mu_p}(s, r_0) = \int_0^\infty e^{-st} \mathcal{F}_{\mu_p}(t, r_0) dt \quad (C1) $$

Consequently, $\tilde{\mathcal{F}}_{\mu_p}(s = 0, r_0)$ is the probability of having reached the edge in the interval $[0, \infty[$. The probability of never having reached the edge at $t$ is then

$$ 1 - \int_0^t \mathcal{F}_{\mu_p}(t, r_0) \quad (C2) $$

If the protein has not touched the edge between $[0, t]$, the probability that it will be phosphorylated between $t$ and $t+dt$ is $\mu_p e^{-\mu_p t} dt$. The probability that it will be phosphorylated for the first time between $t$ and $t+dt$ is therefore

$$ \mu_p e^{-\mu_p t} dt \left(1 - \int_0^t du \mathcal{F}_{\mu_p=0}(u, r_0)\right) \quad (C3) $$

The probability to complete the process at time $t$ is therefore

$$ \mathcal{P}(t, r_0) = \mu_p e^{-\mu_p t} \left[ 1 - \int_0^t du \mathcal{F}_{\mu_p=0}(u, r_0) \right] + \mathcal{F}_{\mu_p}(t, r_0) \quad (C4) $$

where we take $\mu_p = 0$ in the integrant, because if the protein has touched the edge, it cannot have been phosphorylated. And after averaging over $0 \le r_0 \le R$

$$ \mathcal{P}(t) = \mu_p e^{-\mu_p t} \left[ 1 - \int_0^t du \mathcal{F}_{\mu_p=0}(u) \right] + \mathcal{F}_{\mu_p}(t) \quad (C5) $$

To check that this probability is suitably normalized, we note that we have

$$ \tilde{\mathcal{F}}_{\mu_p=0}(s = \mu_p) = \tilde{\mathcal{F}}_{\mu_p}(s = 0) \quad (C6) $$

because the combination $\mu_p + s$ comes into play. We have, indeed

$$ \tilde{\mathcal{F}}_{\mu_p}(s) = g(x) = \frac{2 I_1(x)}{x I_0(x)}, \quad x = R_a \left[ \frac{\mu_p + s}{D} \right]^{1/2} \quad (C7) $$

In (C4), integrating over $t$, we have

$$ \int_0^\infty P(t) dt = \mu_p \int_0^\infty e^{-\mu_p t} dt = 1 \quad (C8) $$

because

$$ \int_0^\infty \mu_p e^{-\mu_p t} \int_0^t du \mathcal{F}_{\mu_p=0}(u) dt = - \tilde{\mathcal{F}}_{\mu_p}(s=0) \quad (C9) $$

The probability of the protein being phosphorylated between $t$ and $t+dt$ is therefore

$$ \mathcal{P}_p(t, l/R) dt = \frac{1}{\mathcal{Z}} \mu_p e^{-\mu_p t} \left[ 1 - \int_0^t du \mathcal{F}_{\mu_p=0}(u) \right] dt \quad (C10) $$

where $\mathcal{Z}$ is a normalization factor.

$$ \mathcal{Z} = 1 - g(1/x) \quad x = l_p/R_a \quad (C11) $$

Its first moment is calculated as

$$ \langle t_{phos} \rangle = \frac{1}{\mathcal{Z}} \int_0^\infty t \mu_p e^{-\mu_p t} \left[ 1 - \int_0^t du \mathcal{F}_{\mu_p=0}(u) \right] dt \quad (C12) $$

This gives for the first moment

$$ \langle t_{phos} \rangle = \frac{1}{\mu_p} \left[ 1 - g(1/x) + \frac{R}{2l_p} g'(1/x) \right] / (1 - g(1/x)) \quad (C13) $$

If we define the diffusion time as

$$ R_a^2 = 8Dt_{diff} \quad \text{or} \quad \mu_p t_{diff} = R_a^2/8l_p^2 \quad (C14) $$

then (C13) gives the solid curve of Fig. 4. As the size of the adhesion disk is finite, the time distribution is cut off and, for proteins that are phosphorylated, the mean time is less than $1/\mu_p$.

Finally, Taylor expanding $g(x)$ for $x \ll 1$ gives $\langle t_{phos} \rangle = 1/\mu_p$ as expected.

# Appendix D: Time spent in the activation disk for a mortal protein

We define a characteristic time $t_{\mu_p}^{in} = 0$ as the mean time that a mortal protein spends in the disk before desorbing from the membrane. From the survival probability calculated above, this reference time can be calculated as follows

$$ 2\pi \int_0^R dr r \int_0^\infty dt \mathcal{P}_1(r,t) \quad (D1) $$

where the Laplace transform of $\mathcal{P}_1(r,t)$ is given in (B10). For $\mu_p = 0$, one has

$$ Z_1 = 1/(qR_a^2) \quad q = \sqrt{\mu_d/D} \quad (D2) $$

Hence

$$ t_{\mu_p=0}^{in} = \frac{1}{\mu_d} \left[ 1 - 2R_a/l_d I_1(R_a/l_d) K_1(R_a) \right] \quad (D3) $$

which allows us to plot the renormalized averaged times for a protein to be phosphorylated on 10 residues in both cases (sequential and random network) as shown in Fig. 10.
## Appendix E: Mean time (model B)

In this appendix, we give the formulas for calculating the average monophosphorylation time in the context of model B. Using (B13) for the Laplace transform, we have

$$ \langle t_{\text{mean}} \rangle = -\mu_p \left[ \frac{1}{\tilde{\mathcal{P}}_1(B)} \frac{\partial \tilde{\mathcal{P}}_1(B)}{\partial s} \right]_{s=0} \quad (\text{E1}) $$

which gives to leading order in the high phosphorylation regime

$$ \langle t_{\text{mean}} \rangle \approx \frac{1}{2\mu_p} \frac{Q(R/l_d)}{K_0(R/l_d)^2} \quad (\text{E2}) $$

with corrections of order $1/\mu_p^{3/2}$ and where

$$ Q(x) = K_0(x)^2 + K_0(x) \left( \frac{2}{x} K_1(x) - K_2(x) \right) + 2K_1(x)^2 \quad (\text{E3}) $$

As expected, $\langle t_{\text{mean}} \rangle \approx 1/\mu_p$ in the small desorption limit, but increases significantly in the other limit.

In the same limit of large phosphorylation rate, one finds for the probability

$$ \mathcal{P}_1(B) \approx 1 - \frac{1}{\mu_p} \left[ 2 \frac{D}{R_d l_d} \frac{K_1(R_d/l_d)}{K_0(R_d/l_d)} + \mu_d \right] \quad (\text{E4}) $$

so that Eqs. (E2) and (E4) give the parametric plot of Fig. 7.

[1] F. Marks, U. Klingmüller, and K. Müller-Decker, Cellular signal processing: An introduction to the molecular mechanisms of signal transduction (Garland Science, New York, 2009).

[2] P. Cohen, Trends Biochem Sci 25, 596 (2000).

[3] M. Köivomügi, M. Ord, A. Iofik, E. Valk, R. Venta, I. Faustova, R. Kivi, E. R. M. Balog, S. M. Rubin, and M. Loog, Nat Struct Mol Biol 20, 1415 (2013).

[4] R. S. Gormal, R. Martinez-Marmol, A. J. Brooks, and F. A. Meunier, Elife 13, e93902 (2024).

[5] A. Remorino, S. De Beco, F. Cayrac, F. Di Federico, G. Cornil-leau, A. Gautreau, M. C. Parrini, J.-B. Masson, M. Dahan, and M. Coppey, Cell Rep 21, 1922 (2017).

[6] Y. Zhou and J. F. Hancock, Biochim Biophys Acta 1853, 841 (2015).

[7] M. F. Garcia-Parajo, A. Cambi, J. A. Torreno-Pina, N. Thompson, and K. Jacobson, J Cell Sci 127, 4995 (2014).

[8] J. Rossy, D. J. Williamson, and K. Gaus, Front Immunol 3, 167 (2012).

[9] B. N. Kholodenko, J. F. Hancock, and W. Kolch, Nat Rev Mol Cell Biol 11, 414 (2010).

[10] D. Bray, Annu Rev Biophys Biomol Struct 27, 59 (1998).

[11] R. Changede, X. Xu, F. Margadant, and M. P. Sheetz, Dev Cell 35, 614 (2015).

[12] M. Bachmann, A. Skripka, K. Weißenbruch, B. Wehrle-Haller, and M. Bastmeyer, J Cell Sci 135, jcs258764 (2022).

[13] H. C. Berg and E. M. Purcell, Biophys J 20, 193 (1977).

[14] D. Shoup and A. Szabo, Biophys J 40, 33 (1982).

[15] S. Redner, A guide to first-passage processes (Cambridge University Press, 2001).

[16] C. Collinet and T. Lecuit, Nat Rev Mol Cell Biol 22, 245 (2021).

[17] D. Kirby, J. Rothschild, M. Smart, and A. Zilman, Phys Rev E 103, 042401 (2021).

[18] N. E. Hynes, P. W. Ingham, W. A. Lim, C. J. Marshall, J. Massagué, and T. Pawson, Nat Rev Mol Cell Biol 14, 393 (2013).

[19] J. E. Toettcher, O. D. Weiner, and W. A. Lim, Cell 155, 1422 (2013).

[20] P. E. Farahani, E. H. Reed, E. J. Underhill, K. Aoki, and J. E. Toettcher, Annu Rev Biomed Eng 23, 61 (2021).

[21] A. Kerjouan, C. Boyault, C. Oddou, E. Hiriart-Bryant, A. Gri-

chine, A. Kraut, M. Pezet, M. Balland, E. Faurobert, I. Bonnet, Y. Coute, B. Fourcade, C. Albiges-Rizo, and O. Destaing, J Cell Sci 134, 10.1242/jcs.254599 (2021).

[22] J. Ju, H. N. Lee, L. Ning, H. Ryu, X. X. Zhou, H. Chun, Y. W. Lee, A. I. Lee-Richerson, C. Jeong, M. Z. Lin, and J. Seong, Cell Reports 40, 10.1016/j.celrep.2022.111080 (2022).

[23] J. de Seze, M. Bongaerts, B. Boulevard, and M. Coppey, Elife 12, RP93180 (2025).

[24] J. E. Toettcher, O. D. Weiner, and W. A. Lim, Cell 155, 1422 (2013).

[25] W. Zhou, W. Li, S. Wang, B. Salovska, Z. Hu, B. Tao, Y. Di, U. Punyamurtula, B. E. Turk, W. C. Sessa, and Y. Liu, Nature Communications 14, 3803 (2023).

[26] J. M. G. Vilar, R. Jansen, and C. Sander, PLoS Comput Biol 2, e3 (2006).

[27] Y. E. Antebi, J. M. Linton, H. Klumpe, B. Bintu, M. Gong, C. Su, R. McCardell, and M. B. Elowitz, Cell 170, 1184 (2017).

[28] L. Saiz and J. M. G. Vilar, Mol Syst Biol 2, 2006.0024 (2006).

[29] C. Albigès-Rizo, O. Destaing, B. Fourcade, E. Planus, and M. R. Block, J Cell Sci 122, 3037 (2009).

[30] D. Shoup and A. Szabo, Journal of Electroanalytical Chemistry and Interfacial Electrochemistry 140, 237 (1982).

[31] P. S. Swain and E. D. Siggia, Biophys J 82, 2928 (2002).

[32] C. Salazar and T. Höfer, The FEBS Journal 276, 3177 (2009), https://febs.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1742-4658.2009.07027.x.

[33] P. François and G. Altan-Bonnet, Journal of Statistical Physics 162, 1130 (2016).

[34] U. Alon, An Introduction to Systems Biology: Design Principles of Biological Circuits, second edition ed. (CRC Press, 2020).

[35] T. W. McKeithan, Proc Natl Acad Sci U S A 92, 5042 (1995).

[36] J. J. Hopfield, Proc Natl Acad Sci U S A 71, 4135 (1974).

[37] J. Ninio, Biochimie 57, 587 (1975).

[38] D. Kirby and A. Zilman, Phys. Rev. E 111, 024408 (2025).

[39] I. V. Gopich and A. Szabo, Protein Sci 25, 244 (2016).

[40] I. V. Gopich and A. Szabo, Biopolymers 115, e23533 (2024).

[41] K. Takahashi, S. Tanase-Nicola, and P. R. ten Wolde, Proc Natl Acad Sci U S A 107, 2473 (2010).

[42] A. Mugler, A. G. Bailey, K. Takahashi, and P. R. ten Wolde,
Biophys J **102**, 1069 (2012).

[43] I. V. Gopich and A. Szabo, Proc Natl Acad Sci U S A **110**, 19784 (2013).

[44] A. Mugler, F. Tostevin, and P. R. ten Wolde, Proc Natl Acad Sci U S A **110**, 5927 (2013).

[45] W. Y. C. Huang, Q. Yan, W.-C. Lin, J. K. Chung, S. D. Hansen, S. M. Christensen, H.-L. Tu, J. Kuriyan, and J. T. Groves, Proc Natl Acad Sci U S A **113**, 8218 (2016).

[46] W. Y. C. Huang, S. Alvarez, Y. Kondo, Y. K. Lee, J. K. Chung, H. Y. M. Lam, K. H. Biswas, J. Kuriyan, and J. T. Groves, Science **363**, 1098 (2019).

[47] L. B. Case, X. Zhang, J. A. Ditlev, and M. K. Rosen, Science **363**, 1093 (2019).

[48] D. Oh, M. Ogiue-Ikeda, J. A. Jadwin, K. Machida, B. J. Mayer, and J. Yu, Proc Natl Acad Sci U S A **109**, 14024 (2012).

[49] W. Bialek, Stability and noise in biochemical switches (2000).

[50] P. B. Warren and P. R. ten Wolde, J Phys Chem B **109**, 6812 (2005).

[51] S. Erlendsson and K. Teilum, Front Mol Biosci **7**, 615565 (2020).

[52] L. Albertazzi, F. J. Martinez-Veracoechea, C. M. A. Leenders, I. K. Voets, D. Frenkel, and E. W. Meijer, Proc Natl Acad Sci U S A **110**, 12203 (2013).

[53] T. Curk, J. Dobnikar, and D. Frenkel, Proc Natl Acad Sci U S A **114**, 7210 (2017).

[54] G. V. Dubacheva, T. Curk, D. Frenkel, and R. P. Richter, J Am Chem Soc **141**, 2577 (2019).

[55] W. A. Lim, Curr Opin Struct Biol **12**, 61 (2002).

[56] X. Xu, T. Masubuchi, Q. Cai, Y. Zhao, and E. Hui, Elife **10**, e74276 (2021).

[57] D. T. Gillespie, Annual Review of Physical Chemistry **58**, 35 (2007).

[58] R. Erban and S. J. Chapman, Phys Biol **6**, 046001 (2009).

[59] O. Rossier and G. Giannone, Exp Cell Res **343**, 28 (2016).

[60] T. Orré, A. Joly, Z. Karatas, B. Kastberger, C. Gabriel, R. T. Böttcher, S. Lévêque-Fort, J.-B. Sibarita, R. Fässler, B. Wehrle-Haller, O. Rossier, and G. Giannone, Nat Commun **12**, 3104 (2021).

[61] B. Meerson and S. Redner, Phys Rev Lett **114**, 198101 (2015).

[62] P. L. Krapivsky and S. Redner, The Journal of Chemical Physics **147**, 214903 (2017), https://doi.org/10.1063/1.5009425.

[63] A. M. Berezhkovskii, L. Dagdug, and S. M. Bezrukov, The Journal of Chemical Physics **150**, 064107 (2019), https://doi.org/10.1063/1.5083808.

[64] B. Fourcade, Physical Review E **101**, 10.1103/Phys-RevE.101.012420 (2020).

[65] J. A. McCammon and S. H. Northrup, Nature **293**, 316 (1981).

[66] R. Zwanzig, Acc. Chem. Res. **23**, 148 (1990).

[67] I. V. Gopich and A. Szabo, J Phys Chem B **120**, 8080 (2016).

[68] I. Bena, International Journal of Modern Physics B **20**, 2825 (2006).

[69] N. R. Smith and O. Farago, Phys. Rev. E **106**, 054118 (2022).

[70] J. Ninio, Proc Natl Acad Sci U S A **84**, 663 (1987).

[71] S. C. Kou, B. J. Cherayil, W. Min, B. P. English, and X. S. Xie, J Phys Chem B **109**, 19068 (2005).

[72] P. Li, S. Banjade, H.-C. Cheng, S. Kim, B. Chen, L. Guo, M. Llaguno, J. V. Hollingsworth, D. S. King, S. F. Banani, P. S. Russo, Q.-X. Jiang, B. T. Nixon, and M. K. Rosen, Nature **483**, 336 (2012).

[73] Y. Sawada, M. Tamada, B. J. Dubin-Thaler, O. Cherniavskaya, R. Sakai, S. Tanaka, and M. P. Sheetz, Cell **127**, 1015 (2006).

[74] J. A. Ubersax and J. E. Ferrell Jr, Nature Reviews Molecular Cell Biology **8**, 530 (2007).

[75] A.-L. Le Roux, M. A. Busquets, F. Sagués, and M. Pons, Colloids Surf B Biointerfaces **138**, 17 (2016).

[76] A.-L. Le Roux, I.-L. Mohammad, B. Mateos, M. Arbesú, M. Gairí, F. A. Khan, J. M. C. Teixeira, and M. Pons, iScience **12**, 194 (2019).

[77] M. Anselmi and J. S. Hub, Sci Rep **10**, 18530 (2020).

[78] R. Erban, S. J. Chapman, and P. K. Maini, A practical guide to stochastic simulations of reaction-diffusion process, arxiv preprint arxiv:0704.1908.

[79] Y. Cao and R. Erban, Bull Math Biol **76**, 3051 (2014).

[80] Y. Cao, P. Linda, and E. Seitaridou, Bulletin of Mathematical Biology **81**, 2819 (2019).

[81] S. S. Andrews, Bioinformatics **33**, 710 (2016).

[82] The dislin library used in this work is available at: http://www.dislin.de/.