Perfect adaptation in eukaryotic gradient sensing using cooperative allosteric binding

TL;DR

This paper models eukaryotic gradient sensing as a cooperative allosteric process, showing cells can adapt to various ligand concentrations and achieve near-optimal chemotaxis, with accuracy influenced by diffusion and reaction rates.

Contribution

It introduces a novel allosteric model for gradient sensing that explains adaptation and accuracy trade-offs in eukaryotic chemotaxis.

Findings

Cells can adapt receptor affinity to achieve 50% receptor binding across concentrations.

Chemotactic accuracy approaches optimal over a broad concentration range.

Diffusion of allosteric compounds significantly impacts sensing accuracy.

Abstract

Eukaryotic cells generally sense chemical gradients using the binding of chemical ligands to membrane receptors. In order to perform chemotaxis effectively in different environments, cells need to adapt to different concentrations. We present a model of gradient sensing where the affinity of receptor-ligand binding is increased when a protein binds to the receptor's cytosolic side. This interior protein (allosteric factor) alters the sensitivity of the cell, allowing the cell to adapt to different ligand concentrations. We propose a reaction scheme where the cell alters the allosteric factor's availability to adapt the average fraction of bound receptors to 1/2. We calculate bounds on the chemotactic accuracy of the cell, and find that the cell can reach near-optimal chemotaxis over a broad range of concentrations. We find that the accuracy of chemotaxis depends strongly on the diffusion of the allosteric compound relative to other reaction rates. From this, we also find a trade-off between adaptation time and gradient sensing accuracy.

Figures

Tables

Table 1. Table 1: Table of simulation parameters. 2 2 2 These parameters are used throughout the paper; any deviation from them is explicitly noted. Description of the parameter and their dimensions are given along with their value in “dimensionless” units—for example, the values of G ~ \tilde{G} , K ~ M \tilde{K}_{M} in Eq. ( 14 ).

Parameter	Description	Unit	Dimensionless value
$g$	Gradient steepness	$1$	$0.05$
$C_0$	Mean concentration of ligands	nM
$N$	Number of receptors	$1$	$4\times {10}^5$
$\Delta \theta$	Angular spacing between simulated G-protein sites	1	$\pi /20$
$K_D$	Dissociation constant with ligand	nM	$100$
$\alpha$	Receptor cooperativity	$1$	$10$
$K_G$	Dissociation constant with allosteric factor	nM	$1$
$K_M$	Michaelis constant	nM	${10}^{-4}\left(=K_M/K_G\right)$
$V_{\text{max}}$	Limiting rate of allosteric protein activation	nM/s	1
$D$	Diffusion coefficient of allosteric proteins	$\text{µ}{m}^2/s$
$\Delta t$	Simulation time step	s	${10}^{-5}\left(=\Delta t\frac{V_{\max }}{K_G}\right)$
$G_{\mathrm{tot}}$	Total allosteric protein concentration	nM	$50\left(=G_{\mathrm{tot}}/K_G\right)$

Equations

$$C(\theta )=C_0\exp \left[\frac{g}{2}\cos (\theta -\varphi )\right],$$

$$\mathcal{L}\left(g,\varphi |\{X_i\}\right)=\prod _{i=1}^N{p}_{\text{bound},i}^{X_i}(1-p_{\text{bound},i}{)}^{1-X_i}.$$

$$\text{Var}({\hat{\theta }}_{\mu })⩾[I(\mathbf{\theta }{)}^{-1}{]}_{\mu \mu },$$

$$I_{\varphi \varphi }=\frac{N{g}^2{C}_0{K}_D}{8(C_0+K_D{)}^2}.$$

$$I_{\varphi \varphi }=\frac{N{\xi }_1^2{g}^2}{2{\xi }_0(1-{\xi }_0)}+\mathcal{O}(g^3).$$

$$p_{\text{bound}}=\frac{C}{C+K_{\text{eff}}},$$

$$K_{\text{eff}}=K_D\frac{G+K_G}{\alpha G+K_G},$$

$$I_{\varphi \varphi }=\frac{N{g}^2{C}_0{K}_{\text{eff}}}{8(C_0+K_{\text{eff}}{)}^2}.$$

$$I_{\varphi \varphi }^{\text{opt}}=\frac{N{g}^2}{8}\{\begin{array}{ll}\frac{C_0{K}_D/\alpha }{(C_0+K_D/\alpha {)}^2},& C_0\le K_D/\alpha ,\\ 1/4,& K_D/\alpha \le C_0\le K_D,\\ \frac{C_0{K}_D}{(C_0+K_D{)}^2},& C_0\ge K_D.\end{array}$$

$$I_{\varphi \varphi }^{\text{2-type}}=\frac{N{g}^2}{16}\left(\frac{C_0{K}_D/\alpha }{(C_0+K_D/\alpha {)}^2}+\frac{C_0{K}_D}{(C_0+K_D{)}^2}\right).$$

$$\frac{\partial G}{\partial t}=V_{\text{max}}\left[-f_b(G,\theta )\frac{G}{G+K_M}+f_u(G,\theta )\frac{G^{\ast }}{G^{\ast }+K_M}\right],$$

$$\frac{\partial G^{\ast }}{\partial t}=V_{\text{max}}\left[f_b(G,\theta )\frac{G}{G+K_M}-f_u(G,\theta )\frac{G^{\ast }}{G^{\ast }+K_M}\right],$$

$$\frac{\partial G}{\partial t}=D{\nabla }^{2}G+\mathcal{R}(G,G^{\ast }),$$

$$\frac{\partial \tilde{G}}{\partial \tilde{t}}=\tilde{D}\frac{{\partial }^2\tilde{G}}{\partial {\theta }^2}-f_b(\tilde{G},\theta )\frac{\tilde{G}}{\tilde{G}+{\tilde{K}}_M}+f_u(\tilde{G},\theta )\frac{\widetilde{G^{\ast }}}{\widetilde{G^{\ast }}+{\tilde{K}}_M},$$

$$I_{\varphi \varphi }\approx \{\begin{array}{ll}\frac{N{g}^2(\alpha -1{)}^2{C}_0^2{K}_D^2{\tilde{D}}^2}{8{\left(\alpha C_0-K_D\right)}^4},& \tilde{D}\ll 1,\\ \frac{N{g}^2}{32}{\left(1-\frac{(\alpha C_0-K_D{)}^2}{2\tilde{D}(\alpha -1)C_0{K}_D}\right)}^2,& \tilde{D}\gg 1.\end{array}$$

$${\tau }_A\equiv \frac{{\int }_{t_0}^{\infty }(t-t_0)\left|I_{\varphi \varphi }(t)-I_{\varphi \varphi }^{\text{ss}}\right|\mathrm{d}t}{{\int }_{t_0}^{\infty }\left|I_{\varphi \varphi }(t)-I_{\varphi \varphi }^{\text{ss}}\right|\mathrm{d}t}.$$

$$\frac{\mathrm{d}{p}^{\circ }}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{p}^{\bullet }}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{p}_G^{\circ }}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{p}_G^{\bullet }}{\mathrm{d}t}$$

$$\left[\begin{array}{cccc}1& 1& 1& 1\\ k_{1}C& -k_2-k_{3}G& 0& k_4/\alpha \\ k_{3}G& 0& -k_4-k_{1}C& k_2/\alpha \\ 0& k_{3}G& k_{1}C& -\frac{k_2+k_4}{\alpha }\end{array}\right]\left[\begin{array}{c}{p}_{}^{\circ }\\ p_{}^{\bullet }\\ p_G^{\circ }\\ p_G^{\bullet }\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right].$$

$$p^{\bullet }$$

$$p_G^{\bullet }$$

$$p_{\mathrm{bound}}=p^{\bullet }+p_G^{\bullet }=\frac{C}{C+K_D\frac{G+K_G}{\alpha G+K_G}}.$$

$$p_{\text{bound}}={\xi }_0+{\xi }_{1}g\cos (\theta -\varphi ),$$

$$\mathcal{L}\left(g,\varphi \mid \{X_i\}\right)=\prod _{i=1}^N{\left({\xi }_0+{\xi }_{1}g\cos ({\theta }_i-\varphi )\right)}^{X_i}\times {\left(1-{\xi }_0-{\xi }_{1}g\cos ({\theta }_i-\varphi )\right)}^{1-X_i},$$

$$\frac{{\partial }^2\ln \mathcal{L}}{\partial {\varphi }^2}=\sum _{i=1}^N\left[X_i(-\frac{g^2{\xi }_1^2{\sin }^2({\theta }_i-\phi )}{({\xi }_0+g{\xi }_1\cos ({\theta }_i-\phi ){)}^2}-\frac{g{\xi }_1\cos ({\theta }_i-\phi )}{{\xi }_0+g{\xi }_1\cos ({\theta }_i-\phi )}\right)$$

$$-(1-X_i)\left(\frac{g{\xi }_1\cos ({\theta }_i-\phi )}{1-{\xi }_0-g{\xi }_1\cos ({\theta }_i-\phi )}+\frac{g^2{\xi }_1^2{\sin }^2({\theta }_i-\phi )}{(1-{\xi }_0-g{\xi }_1\cos ({\theta }_i-\phi ){)}^2})\right].$$

$$I_{\varphi \varphi }$$

$$\frac{\partial \tilde{G}}{\partial \tilde{t}}=\tilde{D}{\partial }_{\theta }^2\tilde{G}-f_b(\tilde{G},\theta )+f_u(\tilde{G},\theta ),$$

$$\tilde{G}={\tilde{G}}_{1/2}(\theta )+ϵ{\tilde{G}}_1(\theta )+\mathcal{O}({ϵ}^2),$$

$$ϵ{\partial }_{\theta }^2{\tilde{G}}_{1/2}-2ϵ{\tilde{G}}_1{\frac{\partial f_b}{\partial \tilde{G}}|}_{\tilde{G}={\tilde{G}}_{1/2}}=0.$$

$$\tilde{G}_{1}={\partial_{\theta}^{2}\tilde{G}_{1/2}}\Big{/}{2\left.\frac{\partial f_{b}}{\partial\tilde{G}}\right|_{\tilde{G}_{1/2}}}.$$

Full text

I Introduction

Eukaryotic cells can sense and follow chemical gradients (a process called “chemotaxis”); this is essential to processes like wound healing, immune response, and cancer metastasis [1, 2, 3]. Eukaryotic chemotaxis has been shown to have remarkable sensitivity to shallow chemical gradients, with some studies finding that cells can sense gradients on the order of 1–2% [4, 5] over a cell width of a few microns. In eukaryotes, chemotaxis occurs by cells recording concentration measurements at different points in space, using membrane receptors to which ligands in the surrounding fluid are bound, creating downstream signals of cell polarity and motility  [6]. Since the biochemistry of ligand-receptor binding is subject to noise via thermal processes, the resulting signal is stochastic in nature [7, 5, 8]. In many circumstances, if the stochasticity arising from ligand-receptor binding limits the ability of a cell to sense a gradient, gradient sensing accuracy is optimal in a narrow concentration range near the ligand-receptor dissociation constant [9, 5, 10, 11, 12]. This narrow range may be a problem for cells whose natural environments are challenged with a broad range of chemoattractant concentrations [13]. Is there a mechanism where cells can keep their sensing accuracy over wide concentration ranges? One possibility, motivated by the observation of multiple receptor types for a single chemoattractant [14, 15], is that cells may hedge their bets by expressing multiple receptor types with different dissociation constants [16] but in doing so sacrifice peak sensing accuracy. Another possibility, which we study here, is an adaptive modification to receptors to maintain chemotactic accuracy. Receptor-level adaptation is well-known in bacterial systems but not believed to be present in eukaryotes [17] – eukaryotic chemotaxis has adaptation mechanisms downstream of receptor binding [18]. In this paper, we propose a model that allows eukaryotic cells to effectively adapt ligand-receptor affinity through the regulation of a protein that allosterically changes the dissociation constant of the receptor. This could include a G protein, arrestin, or any allosteric modulator well described by a ternary complex model [19, 20, 21, 22]. If the allosteric protein’s availability is regulated by the local fraction of receptors bound, the cell can control dissociation constants to ensure the fraction bound is close to 1/2, where chemotaxis is most accurate. We find that near-perfect adaptation is possible under certain constraints on the parameters involved (Secs. II.2–II.4). We claim that the diffusion coefficient of this allosteric protein relative to its activation rate, is what decides whether such adaptation is advantageous for the cell, which we evaluate by looking at the gain in accuracy over simple single-receptor-type chemotaxis or bet-hedging. However, if the activation rate is slow, allowing near-perfect adaptation, this comes with an unavoidable tradeoff of a long time required to adapt to new concentrations. This may explain why, even though eukaryotic cells have all the necessary requirements to perfectly adapt at the receptor level, no evidence of this has been seen.

II Model and results

We review how a cell’s ability to sense a chemical gradient is related to the stochasticity of its receptor-ligand interactions in Sec. II.1. We then use these results to understand how gradient sensing is modified depending on the concentration of an allosteric protein in a model where it can enhance ligand-receptor binding in Sec. II.2. We then propose in Sec. II.3 a scheme where cells can adapt their concentration of active allosteric proteins to maximize the cell’s ability to sense a chemical gradient’s orientation. In Sec. II.4 we show this scheme requires sufficiently fast diffusion of allosteric protein to lead to perfect adaptation, and show in Sec. II.5 that this constraint implies a tradeoff between accuracy of sensing and time required to adapt to new concentration levels.

II.1 Gradient sensing accuracy depends on ligand binding probabilities

We model cells in two dimensions as circular with receptors on the boundary. We assume cells are in an exponential concentration gradient,

[TABLE]

where $C_{0}$ is the mean concentration and $g$ the percentage change across the cell, $R$ is the cell radius, $\theta$ is the angular coordinate of a receptor and $\phi$ is the gradient direction. We will later assume shallow gradients $g\ll 1$ to make certain analytical calculations tractable. The state of the $i$th receptor can be modeled as a Bernoulli random variable $X_{i}$, where $X_{i}=1$ indicates the receptor is bound and $X_{i}=0$ indicates the receptor is unbound. If there are $N$ receptors on the cell membrane, the associated likelihood function is:

[TABLE]

We can then extract the Cramér-Rao bound for the variance of the unbiased estimators $\hat{\bm{\theta}}\equiv(\hat{g},\hat{\phi})$ of maximum likelihood estimation (MLE) [10, 23]:

[TABLE]

where $I(\bm{\theta})^{-1}$ denotes the inverse of the Fisher information matrix $I_{\mu\nu}=-\langle\partial^{2}\ln\mathcal{L}/\partial\theta_{\mu}\partial\theta_{\nu}\rangle$, where the indices $\mu,\nu$ run over the gradient magnitude $g$ and orientation $\phi$. This gives a lower bound for the variance of any unbiased estimator – this reflects the best estimate that a cell can make given the information it has. Hu et al. [10] calculated these bounds for a 2D circular cell with the concentration profile of Eq. (1), and assuming that there are simple ligand-receptor kinetics. In this case, the probability of being bound is $p_{\textrm{bound},i}={C(\theta_{i})}/{\left[C(\theta_{i})+K_{D}\right]}$, where $K_{D}$ is the ligand-receptor dissociation constant. The $(\phi,\phi)$ term in the Fisher information matrix is then [10]

[TABLE]

We think of $I_{\phi\phi}$ as characterizing the minimal error possible in a cell’s estimation of the direction of the gradient – the Cramér-Rao bound for the gradient direction is $\textrm{Var}(\hat{\phi})\geqslant I^{-1}_{\phi\phi}$. Sweeping over the mean concentration $C_{0}$, the Fisher Information peaks at $C_{0}=K_{D}$, and falls off monotonically on either side of the peak.

We can generalize Eq. (4) to an arbitrary shallow-gradient concentration profile by expanding $p_{\textrm{bound}}$ in terms of $g$. Given a bound probability $p_{\textrm{bound}}(\theta)=\xi_{0}+\xi_{1}g\cos(\theta-\phi)+\mathcal{O}(g^{2})$, we find $I_{\phi\phi}$ as

[TABLE]

The details of this calculation can be found in Appendix B. This proves to be a useful result as it allows for a quick evaluation of gradient sensing accuracy given the spatial profile of receptor binding.

II.2 Cooperative ternary interactions allow the cell to tune its region of maximal accuracy

We model the allosteric protein-receptor interaction using a ternary complex model [21] involving receptor, ligand and allosteric protein. Each receptor can bind one ligand molecule on the extracellular side, and one protein molecule on the cytosolic side (Fig. 1). The dissociation constants for these processes are $K_{D}$ for ligand binding to the bare receptor and $K_{G}$ for the allosteric protein binding to the bare receptor. When this protein binds to the receptor, this decreases the ligand-receptor dissociation constant from $K_{D}$ to $K_{D}/\alpha$. Unless otherwise noted, we assume $\alpha>1$, i.e., cooperative binding – the ligand-receptor interaction is stronger if the G protein is bound – and take $\alpha$ to be of order 10 [21]. From detailed balance, the dissociation constant associated with the allosteric protein binding must be $K_{G}/\alpha$ when a ligand is bound [24]. In this model, the total probability that a receptor is bound to a ligand molecule can be written as

[TABLE]

with an effective dissociation $K_{\text{eff}}$ given by

[TABLE]

where $G$ is the allosteric protein concentration (Appendix A).

Eq. (7) shows that the effective ligand-receptor interaction can have a dissociation ranging from $K_{D}/\alpha$ (in the limit $G\gg K_{G}$) to $K_{D}$ (in the limit $G\ll K_{G}$). If all receptors on the membrane have bound probability of Eq. (6), we can compute the Fisher information from Eq. (5), finding:

[TABLE]

We find that using the ternary complex model, we get a Fisher information that is exactly that of Eq. (4) but with the effective dissociation constant $K_{\textrm{eff}}$. The Fisher information $I_{\phi\phi}$ is thus maximized when $C_{0}=K_{\text{eff}}$. This point of maximal accuracy is when the probability of binding is 1/2, by Eq. (6). The presence of the allosteric protein allows the cell to tune its region of peak efficiency to any value between $K_{D}/\alpha$ and $K_{D}$ by changing $G$. We show that this can lead to perfect adaptation of sensing accuracy below in Sec. II.3.

If $K_{\textrm{eff}}$ takes on its optimal possible value, i.e. $K_{\textrm{eff}}=C_{0}$ when $K_{D}/\alpha\leq C_{0}\leq K_{D}$, we find

[TABLE]

We plot this optimal Fisher information in Fig. 2. We see that the Fisher information is near its maximum value for a broad range of concentrations.

We compare the optimal adaptation result of Eq. (9) to two possible alternatives. First, a cell with a simple single receptor with dissociation constant $K_{D}$, leads to a clear peak in information near $K_{D}$. Secondly, if instead the cell has $N/2$ receptors with dissociation constant $K_{D}$ and $N/2$ with $K_{D}/\alpha$, this hedges the cell’s bets [16]—this allows cells to chemotax accurately at a broader range of concentrations but at a cost of lowering the maximum Fisher information. For example, a cell with two receptor types with dissociation constants $K_{D}$ and $K_{D}/\alpha$ has Fisher jnformation

[TABLE]

If a cell can actually reach the optimal adaptation of Eq. (9), it will have more information than either of these two alternatives.

There is an implicit assumption behind Eq. (8). In calculating this, we assumed that the cell only “knows” whether or not a receptor is bound, i.e. that from the perspective of the cell, it cannot distinguish between the right two states in Fig. 1). This is naturally an oversimplification of GPCR signaling. For instance, if our allosteric protein were a G protein, in the canonical model of G protein-based signaling, activated G proteins must unbind from the receptor to lead to downstream signaling, something we do not take into account here. Our assumption here is made in large part so that our model will limit back to [9, 10]. Different choices for which states are active/signaling (i.e. which states the cell “knows” about) and which states the cell can distinguish (i.e., whether the likelihood includes a binary $X_{i}$ or a ternary one) can be made and will lead to different results.

II.3 Feedback based on bound ligand fraction leads to perfect adaptation

For the cell to maximize its accuracy, it must change $G$ to tune $K_{\textrm{eff}}$ to the typical concentration of ligand around it, $C_{0}$ – or to regulate the probability of binding to be 1/2. We propose a model where the allosteric protein switches between one state that is capable of binding to the receptor ($G$) and one that cannot ($G^{*}$) (Fig. 3). We assume that the reaction $G\to G^{*}$ is catalyzed either by the bound receptor or by an enzyme downstream of the bound receptor, so the rate is proportional to the bound receptor fraction, and similarly the reverse process catalyzed by unbound receptor is proportional to the unbound receptor fraction. If these rates are of the Michaelis-Menten form [25, 26], then we can write an equation for the concentration of $G$ at a particular angular location $\theta$ along the cell:

[TABLE]

where $K_{M}$ is a Michaelis constant, $G$ and $G^{*}$ are the active and inactive allosteric protein concentrations, $f_{b}$ and $f_{u}=1-f_{b}$ are the local fraction of bound and unbound receptors at $\theta$. We note that, in principle, the fraction bound is a stochastic variable, reflecting ligand-receptor binding fluctuations. We use the $f_{b}$ language to emphasize this potential complication. However, we neglect this stochasticity, assuming the fraction bound is just $f_{b}(G,\theta)\approx\langle f_{b}(G,\theta)\rangle=p_{\textrm{bound}}(G,C(\theta))$. Essentially, when we write this equation, we are writing an equation for a region of angles near $\theta$, and assuming $G$ is well-mixed over this region. Using the deterministic model seems reasonable if the number of receptors in this segment of the cell membrane is large, or the dynamics of $G$ is slow with respect to the kinetics of binding and unbinding, leading to effective time-averaging. Our deterministic model might not be appropriate in the limit of small diffusion coefficients, when the angular range where $G$ can be treated as well-mixed will be small. However, since we find our adaptation mechanism fails in the limit of small $D$ even in the deterministic case, we do not address this further. Because we assume only the active form of this protein can bind to receptors, $f_{b}$ depends on $G$ through Eq. (6). $V_{\text{max}}$ is the maximum rate at which the allosteric protein can be activated/deactivated. The total concentration of these proteins $G_{\textrm{tot}}=G+G^{*}$ is fixed, and so we have a corresponding equation for $G^{*}$:

[TABLE]

In the limit where $G,G^{*}\gg K_{M}$, Eq. (11) reduces to $\frac{\partial G}{\partial t}=V_{\text{max}}\left[-f_{b}+f_{u}\right]$ – so the system will reach a steady state where $f_{u}=f_{b}=1/2$ – leading to perfect adaptation [27]. However, because $K_{\textrm{eff}}$ can only range from $K_{D}/\alpha$ to $K_{D}$, the adaptation will not be able to reach this steady state when the concentration $C$ is outside this range. We also see that perfect adaptation requires $K_{M}$ to be small. We show in Fig. 4(a) a plot of fraction bound [Eq. (6)] as a function of ligand concentration, where we find $G$ [Fig. 4(b)] by simulating Eq. (11) until reaching steady state at each concentration (see Appendix 40). We see that in the range from $K_{D}/\alpha$ to $K_{D}$ that the fraction bound adapts to $\sim 1/2$, and this adaptation is more accurate for smaller Michaelis constant. Outside of the range $K_{D}/\alpha<C<K_{D}$, the concentration of active protein $G$ cannot go above the total concentration $G_{\textrm{tot}}$ or below zero [Fig. 4(b)] and adaptation fails. In the large and small-concentration limits, the binding probability behaves simply as if the receptors had dissociation constant $K_{D}/\alpha$ for $C<K_{D}/\alpha$, or $K_{D}$, for $C>K_{D}$.

We show how the allosteric protein concentration and binding probability adapt to a step change in concentration in Fig. 4(c). The bound probability takes a sudden dip as concentration decreases, but as the concentration of the allosteric protein increases (therefore decreasing $K_{\textrm{eff}}$), the probability of a ligand being bound evolves to near 1/2 again. We note that adaptation of the fraction bound here does not quite reach 1/2, due to the finite value of the Michaelis constant $K_{M}$ (e.g. as seen in [27] for generic negative feedback adaptation).

II.4 Combining adaptation and protein diffusion

In the previous section, we implicitly assumed that the allosteric protein $G$ was localized at a particular angular location $\theta$. Though this scheme allows for perfect adaptation of $p_{\textrm{bound}}$, it is not ideal for gradient sensing, because each receptor will then adapt to the local concentration, resulting in each receptor being bound with probability $\approx 1/2$ where $K_{\text{eff}}(\theta)\approx C(\theta)$. In that case, there is no gradient of bound receptors across the cell (though there will be a gradient of allosteric proteins). To maximize the change in number of bound receptors across the cell, $K_{\text{eff}}$ should be constant across the cell – essentially adapting to the concentration $C_{0}$ at the cell midpoint, instead of the local concentration. This can happen if the allosteric protein can diffuse on the cell membrane. We therefore incorporate a diffusion term into Eq. (11):

[TABLE]

where the reaction term $\mathcal{R}(G,G^{\ast})$ is just the right-hand side of Eq. (11). $G^{\ast}$ then obeys the complementary equation $\partial G^{\ast}/\partial t=D\nabla^{2}G^{\ast}-\mathcal{R}(G,G^{\ast})$, so the total allosteric protein (active+inactive) concentration ($\oint(G+G^{\ast})\mathrm{d}\mathbf{x}$) is conserved. Note that on the 1D surface of a 2D circular cell, $\nabla\equiv(1/R)\partial_{\theta}$. The adaptive process now has natural scales of distance ($R$) and time ($K_{G}/V_{\text{max}}$), so we rescale Eq. (13) to make it dimensionless:

[TABLE]

where $\tilde{G}=G/K_{G}$, $\tilde{K}_{M}=K_{M}/K_{G}$, $\tilde{t}=tV_{\text{max}}/K_{G}$ and $\tilde{D}=DK_{G}/V_{\text{max}}R^{2}$. When diffusion is slow relative to the switching between states of the allosteric protein, $\tilde{D}\ll 1$, we expect that the bound fraction will be approximately 1/2 everywhere –i.e. the cell has almost no information about the gradient direction. As we increase $\tilde{D}$, we expect the cell’s gradient sensing to approach the optimal case described in Sec. II.2 – in this limit, it is equivalent to all the receptors adapting $K_{\text{eff}}$ to be the mean concentration $C_{0}$.

We solve Eq. (14) numerically, using forward time-centered space (FTCS) method with receptors on the surface of such a cell with radius $R$, using a time step $\Delta\tilde{t}=0.001$ and $\tilde{K}_{M}=0.0001$ in the units described earlier. We simulate Eq. (14) until the equations reach their steady state value $G(\theta)$. Using these values of $G(\theta)$ then gives us a prediction for how the bound probability $p_{\textrm{bound}}$ depends on angle $\theta$. We fit this simulation result to the form $p_{\textrm{bound}}(\theta)=\xi_{0}+\xi_{1}g\cos(\theta-\phi)$ using linear regression, which allows us to predict the Fisher information from Eq. (5). We can obtain steady-state Fisher information matrix elements for varying dimensionless diffusion coefficients $\tilde{D}$ (Fig. 5). At large $\tilde{D}$, we see that cells with the adaptation process nearly reach the optimal possible Fisher information [Fig. 5(a)], as we expect. However, we do note that even at $\tilde{D}=50$ there is a slight decrease in $I_{\phi\phi}$ as $C_{0}$ increases from $K_{D}/\alpha$ to $K_{D}$. By contrast, for low diffusion or fast reactions $\tilde{D}=0.1$, the Fisher information can drop orders of magnitude below the optimal value, far worse than the alternatives shown in Fig. 2. (We note that comparisons between the lowest values of the Fisher information should not be taken too seriously – at very small Fisher information, the Cramer-Rao bound will be misleading because the variance of the angle $\phi$ is bounded because $\phi$ must be between [math] and $2\pi$ [28, 29].)

In addition to our numerical solutions to Eq. (14), we use a perturbative approach to solve for the steady-state allosteric protein concentration. We perturb about the $\tilde{D}=0$ case to order $\tilde{D}$ and around $1/\tilde{D}=0$ to order $1/\tilde{D}$. The full solution is somewhat lengthy and detailed in Appendix C. We obtain the following results for the perturbed Fisher information in the range $C_{0}\in(K_{D}/\alpha,K_{D}]$:

[TABLE]

As we expect, in the limit $\tilde{D}\to 0$, the Fisher information vanishes – in this limit, the receptor binding is uniform across the cell. In the limit $\tilde{D}\gg 1$, we see that the Fisher information converges to the optimal value of $Ng^{2}/32$. For both asymptotic limits, we note that at finite values of $\tilde{D}$, the Fisher information decreases as $C_{0}$ increases from $K_{D}/\alpha$ to $K_{D}$. We believe this happens because when $G$ is high (smaller $C_{0}$), small variations in $G$ do not result in large changes in $K_{\mathrm{eff}}$, and in our model, a flatter profile of $K_{\text{eff}}(\theta)$ over the cell around a mean value of $K_{\mathrm{eff}}=C_{0}$ gives a higher $I_{\phi\phi}$. Fig. 5 shows that the numerical solution agrees with the perturbation results in both limits.

Under what circumstances would our adaptation model be more accurate than a cell using a single receptor type [9], or using two receptor types [16], as we showed in Fig. 2? We plot the difference in Fisher information between the adaptation model and these baseline models $\Delta I_{\phi\phi}=I_{\phi\phi}^{\textrm{adapt}}-I_{\phi\phi}^{\textrm{baseline}}$ in Fig. 6. We see immediately that $\tilde{D}$ must be on the order of 10–100 in order for the adaptation model to outperform both models. Improvements over the two-receptor type model are less dramatic, but the difference in Fisher information $\Delta I_{\phi\phi}$ is positive over a larger region in the parameter space.

We expect that the cooperativity $\alpha$ affects the accuracy of adaptation by determining the degree to which the allosteric protein affects ligand binding – a naive guess would be that a larger cooperativity $\alpha$ would allow adaptation over a larger range of concentrations. However, this is only true in the limit $\tilde{D}\to\infty$. In Fig. 7 we scan over values of $\alpha$ and the mean concentration $C_{0}$ at two different values of diffusion $\tilde{D}=10$ and $\tilde{D}=100$. We do see that at this highest value of $\tilde{D}$ that increasing $\alpha$ increases the range of concentrations $C_{0}$ where near-optimal sensing is possible. However, at lower $\tilde{D}=10$, increasing $\alpha$ just shifts the region of accurate sensing around $K_{D}/\alpha$.

As discussed in Sec. II.3, increasing the Michaelis constant should affect gradient sensing accuracy at higher concentrations near $C_{0}=K_{D}$. This can be seen in Fig. 8. When $K_{M}$ is sufficiently small, we observe a rapid increase in gradient sensing accuracy near $K_{D}$, after which the Fisher information seems to fall of like in the single-type model [Eq. (4)]. For larger $K_{M}$, however, the sharp jump is smoothed out. This follows from the fact that smaller $K_{M}$ makes it easier to maintain the limit $G,G^{*}\gg K_{M}$ where near-perfect adaptation is possible (Sec. II.3).

II.5 Trade-offs in sensing: accuracy vs time to adapt

Our results show that the adaptation model is more successful when $\tilde{D}=DK_{G}/V_{\text{max}}R^{2}$ is large. This can occur either by increasing the diffusion coefficient of the allosteric protein or by decreasing the production rate $V_{\textrm{max}}/K_{G}$. However, because the diffusion coefficient of the protein cannot be increased arbitrarily, we think it is more plausible that $\tilde{D}$ can be increased by decreasing $V_{\textrm{max}}$. However, when the cell decreases $V_{\textrm{max}}$, this means that the conversion between $G$ and $G^{*}$ is slower – so the time required for the cell to adapt to a new concentration will also slow.

In Fig. 9, we simulate gradient sensing and adaptation for a broad range of values of $\tilde{D}$. We treat the physical diffusion coefficient as fixed to $D=0.1\ \textrm{\textmu}\mathrm{m^{2}/s}$ and cell size fixed to $R=5\ \textrm{\textmu}\textrm{m}$ Given these values and a value of $\tilde{D}=DK_{G}/V_{\text{max}}R^{2}$, this lets us set the characteristic timescale $K_{G}/V_{\textrm{max}}$ in real units. We define an adaptation time $\tau_{A}$ as the characteristic timescale for the change in the Fisher information in response to a jump in the concentration. Assuming the jump occurs at time ${t}_{0}$ and the Fisher information relaxes to a steady-state value $I_{\phi\phi}^{\textrm{ss}}=\lim_{t\to\infty}I_{\phi\phi}({t})$, we define the adaptation time as:

[TABLE]

In Fig. 9, we compute this adaptation time for a broad range of $\tilde{D}$, and plot the steady-state Fisher information as a function of this adaptation time. We see generally that higher Fisher information requires longer adaptation times. We also compare the adaptation model to gradient sensing accuracy with fixed receptor types (the blue dash-dot line indicates a 50-50 mix of $K_{D}$ and $K_{D}/\alpha$, the red dashed line a single receptor type of $K_{D}$). The adaptation time required for a particular Fisher information will depend on the final concentration $C_{0}$ (Fig. 9a-c), as we would expect based on the strong dependence of accuracy on $C_{0}$ in Fig. 5. We find that the adaptation time required for the adaptation model to beat the two-receptor-type model is generally of the order of a few hundred seconds. We note that interpreting the adaptation time may not be obvious because the Fisher information may not monotonically vary in adaptation (Fig. 9d). Adaptation times around $\tau_{A}\sim 1000$ seconds might still be useful for slow cells like fibroblasts  [30], though these cells are also larger in size and would have correspondingly larger adaptation times. We show adaptation to relatively modest concentration changes in Fig. 9; we also test $\tau_{A}$ vs accuracy for a different concentration jump size to the same final concentrations (Fig. S1), where we see that the size of the jump in concentration does not significantly affect the adaptation time.

III Discussion

In this work, we have proposed a model of eukaryotic chemotaxis with near-perfect adaptation of receptor properties, allowing cells to sense gradients accurately over a wider range of average concentrations. We found that this approach’s success depends on the relative rates of diffusion and modification of the allosteric protein, and therefore defined the dimensionless diffusion coefficient $\tilde{D}=DK_{G}/V_{\text{max}}R^{2}$. Using a first-order perturbative approach, we derived analytical solutions for the Fisher information in the low and high $\tilde{D}$ limits. We evaluated the effectiveness of the model relative to other previously proposed models of chemotaxis with one/multiple receptor types. We find that adaptation works well for high $\tilde{D}$– but this comes inevitably with a long adaptation time to a change in chemoattractant concentration.

Our model is intentionally simple, and some of the assumptions made are relatively optimistic for the cell’s ability to perform adaptation, and real cells may face significant additional challenges. For instance, we assume in our reaction kinetics [Eq. (11)], that the inactivation rate does not depend on whether the allosteric protein is bound to the receptor. In thermal equilibrium, the rate of transition from $G\to G^{*}$ and its reverse are related to the change in free energy through detailed balance, and the free energy of the state will depend on whether G is bound. Our model is out of equilibrium, and adaptation therefore necessitates an energy cost for the cell, as has been previously proposed [31, 17]. As with the larger program of understanding energy costs for sensing [32], these may further disfavor the adaptation model.

Is adaptation plausible over biologically relevant timescales? For adaptation to have worthwhile results, the cell would need to have $\tilde{D}$ of order $10$, or an adaptation time of on the order of 1000 seconds if $D=0.1\,\textrm{\textmu}\mathrm{m^{2}/s}$ (Fig. 9). Is $\tilde{D}\approx 10$ reasonable? Diffusion coefficients of G-proteins are roughly of order of magnitude $0.1\,\textrm{\textmu}\mathrm{m^{2}/s}$ [33], so the timescale of 1000 seconds is consistent with G proteins or other membrane-attached proteins as the allosteric protein. However, these values imply a $K_{G}/V_{\textrm{max}}$ of $\sim 1000$ s. This is much slower than typical timescales found in phosphorylation and methylation processes in cells [34], but may nonetheless be plausible as such slow timescales have been found in some enzyme activity [35]. Could cells function effectively if this adaptation time were so slow? We compare the adaptation timescale to how fast the cell’s concentration environment changes. If a cell is traveling at a speed of $5\,\textrm{\textmu}\mathrm{m/min}$ in an environment $C(x)=C_{0}e^{gx/L}$ and a $5\%$ gradient across the cell length ($g=0.05$), then a cell must travel $\ln 2/0.05\approx 14$ cell lengths to double the concentration $C_{0}$. This would take 28 minutes if the cell length is 10 $\textrm{\textmu}\mathrm{m}$. This suggests that – in shallow gradients at least – cells may see concentration changes that can be adapted to at relevant timescales. However, if there are rapidly-moving sources, concentration dynamics can be much faster and more complicated [13], and we expect receptor adaptation at the timescale of hundreds or thousands of seconds to not be particularly useful.

The adaptation mechanism we propose here may be more broadly relevant if the allosteric factor has a larger diffusion coefficient, e.g. if it is a primarily cytosolic protein that only occasionally binds to the membrane. (In this case, we would have to more carefully treat a difference in diffusion coefficient between receptor-bound allosteric proteins, which will have low diffusion coefficients, and cytosolic. Within our current model the bound and unbound forms are lumped with a single effective diffusion coefficient.) More speculatively, some GPCRs respond to mechanical cues including membrane tension[36], allowing the membrane itself to act as an “allosteric” factor [37, 38, 39, 40]. Tension could then serve as a global factor [41] in place of an allosteric protein. However, we do note that tension propagation is not straightforward – depending on cell type and on whether membrane or cortex tension is measured, very different extents of tension correlation and propagation can be found [42, 43, 41, 44].

A recent paper also proposed a related idea of using G-protein activation to improve gradient sensing, albeit with different underlying mechanisms. Ghose et al. [45] propose a model of ratiometric sensing where G-proteins serve the role of enhancing signal as well as temporally storing information on ligand from membrane receptors. We do not study the effect of time averaging here, instead focusing on the instantaneous state of whether receptors are bound or not. This is because time averaging in systems where receptor kinetics vary from receptor to receptor may be complicated and a naive average from integrating signals, as we would get from studying the signal-to-noise ratio of $G(t)$, can be very suboptimal [16]. We defer consideration of this point to future work. Surprisingly, given how different our mechanism is from that of the G protein model of [45], Ghose et al. also find a tradeoff between adaptation time and accuracy, with adaptation times on the order of 50-100 seconds. We suspect the tradeoff between ability to adapt to new environments and accuracy in an individual environment is a universal tradeoff, potentially reflecting the time required to communicate information across the cell [46].

Acknowledgements.

The authors acknowledge support from NIH Grant No. R35GM142847 and NSF DMR 1945141. We thank Grace Luettgen and Emiliano Perez Ipiña for useful discussions and a close reading of the draft.

Appendix A Solving the ternary complex model

We write the master equation corresponding to Figure 1, using on/off rates that correspond to the given dissociation constants, as shown in Fig. 10. Without loss of generality, when the dissociation constant is reduced by a factor of $\alpha$, we choose the corresponding off rate to be reduced by a factor of $\alpha$. For the four states, we find the master equation:

[TABLE]

where $C$ and $G$ denote ligand and allosteric protein concentrations, the superscripts $\circ/\bullet$ denote ligand unbound/bound states, and the subscript $G$ indicates that a allosteric protein molecule is bound to the receptor’s cytosolic side. We can take Eqs. (17b)–(17d) in steady state, together with the constraint $p^{\circ}+p^{\bullet}+p_{G}^{\circ}+p_{G}^{\bullet}=1$:

[TABLE]

To be consistent with the notation used in Fig. 1, we define $K_{D}\equiv k_{2}/k_{1}$ and $K_{G}\equiv k_{4}/k_{3}$. Solving for $p^{\bullet}$ and $p_{G}^{\bullet}$ in this steady state, we obtain:

[TABLE]

Then, the total ligand-bound probability should be

[TABLE]

Appendix B Calculating the Fisher information from the bound probability

In this section, we generalize the results of [9, 10] for the Fisher information of a gradient-sensing cell beyond simple ligand-receptor kinetics to any slowly-varying bound probability. If we have bound probability has the form

[TABLE]

where $\xi_{1}g$ is small compared to $\xi_{0}$, we can, assuming Bernoulli statistics as in the main text and [9, 10], form a likelihood function $\mathcal{L}$, finding

[TABLE]

Taking the expected value, which replaces the Bernoulli variables $X_{i}$ with their associated probabilities, we can get the $(\phi,\phi)$ element of the Fisher information matrix to second order in $g$:

[TABLE]

where in the second equality we have approximated the discrete sum as an integral over the cell boundary in the limit of large receptor number $N$.

In fact, we can get the form of this Fisher information up to a geometric prefactor from a simple signal-to-noise calculation akin to [47, 48, 12]. If we simplify our geometry, treating the cell as a line, and have $N/2$ receptors on the front half of the cell and $N/2$ on the back half, the average difference in receptors bound between the back and the front is $\Delta n_{b}=\frac{N}{2}p_{\textrm{bound}}(\textrm{front})-\frac{N}{2}p_{\textrm{bound}}(\textrm{back})=N\xi_{1}g$. If the receptor binding events are independent, then the variance in $\Delta n_{b}$ is $\textrm{Var}(\Delta n_{b})=\textrm{Var}(n_{b}(\textrm{front}))+\textrm{Var}(n_{b}(\textrm{back}))$. These variances are – applying the central limit theorem to the N/2 Bernoulli events of a ligand binding – $\textrm{Var}(\Delta n_{b})=\frac{N}{2}p_{\textrm{bound}}(\textrm{front})(1-p_{\textrm{bound}}(\textrm{front}))+\frac{N}{2}p_{\textrm{bound}}(\textrm{back})(1-p_{\textrm{bound}}(\textrm{back}))$. To the lowest order in $g$, we can neglect the difference between the front and the back, and the probability of being bound is just $\xi_{0}$, so $\textrm{Var}(\Delta n_{b})=N\xi_{0}(1-\xi_{0})$. We thus find that a signal-to-noise ratio $\textrm{SNR}=\Delta n_{b}^{2}/\textrm{Var}(\Delta n_{b})$ is just $\textrm{SNR}=\frac{Ng^{2}\xi_{1}^{2}}{\xi_{0}(1-\xi_{0})}$. We can see that the Fisher information $I_{\phi\phi}$ is proportional to the signal-to-noise ratio.

Appendix C First-order perturbation to allosteric protein concentration

We present perturbative results for both the low- and high-diffusion regimes, which are first-order perturbations for $\tilde{D}\ll 1$ and $\tilde{D}\gg 1$, respectively.

C.1 Low-diffusion regime

If we assume that $G,G^{*}\gg K_{M}$, then the Michaelis terms in Eq. (14) become zero-order, and we find

[TABLE]

where $f_{b}$ is the local bound fraction of receptors. We note that this approximation only holds in the range $K_{D}/\alpha\leq C_{0}\leq K_{D}$, as outside that range, we cannot self-consistently have both $G$ and $G^{*}$ large relative to $K_{M}$. In the limit of $\tilde{D}=0$, the bound fraction will adapt to exactly 1/2 everywhere. We can consider a first-order perturbation from the zero-diffusion case $\tilde{D}=0$:

[TABLE]

where $\epsilon\equiv\tilde{D}$, and $\tilde{G}_{1/2}(\theta)$ is such that $f_{b}\bm{(}\tilde{G}_{1/2}(\theta),\theta\bm{)}=1/2$, which is what we see when $\tilde{D}=0$. The term $f_{b}(\tilde{G},\theta)$ then becomes $f_{b}\approx 1/2+\epsilon\tilde{G}_{1}\left.({\partial f_{b}}/{\partial\tilde{G}})\right|_{\tilde{G}=\tilde{G}_{1/2}}$ and similarly $f_{u}(\tilde{G},\theta)\approx 1/2-\epsilon\tilde{G}_{1}\left.({\partial f_{b}}/{\partial\tilde{G}})\right|_{\tilde{G}=\tilde{G}_{1/2}}$. In steady state ($\partial\tilde{G}/\partial\tilde{t}=0$), to first order in $\epsilon$, we have from Eq. 23,

[TABLE]

from which we obtain

[TABLE]

We can compute both these terms. We can find

[TABLE]

from solving $C(\theta)/(C(\theta)+K_{\textrm{eff}})=1/2$. Taking only terms up to first order in $g$ in this formula, we get

[TABLE]

The local average fraction of bound receptors is $f_{b}={C}/{(C+K_{\textrm{eff}})}$, where the dimensionless $K_{\textrm{eff}}=K_{D}(\tilde{G}+1)/(\alpha\tilde{G}+1)$. Evaluating the other term in Eq. (26):

[TABLE]

we obtain

[TABLE]

where $C$ is given in Eq. (1).

Finally, we have, to first order in $g$,

[TABLE]

Notice the $\cos(\theta-\phi)$ term. If we look at $\tilde{G}_{1/2}(\theta)$, it seems to follow the form $\Gamma_{0}+\Gamma_{1}g\cos(\theta-\phi)$ (at least, to first order in $g$). We can substitute the first-order expansion of $\tilde{G}_{1/2}$ in $g$, then, into the final expression for $G$:

[TABLE]

where $\tilde{G}_{\textrm{opt}}={(K_{D}-C_{0})}/{(\alpha C_{0}-K_{D})}$ is the (uniform) allosteric protein concentration that maximizes the Fisher information when $K_{\textrm{eff}}$ is constant across the cell. Then, plugging this into Eq. (6), we get

[TABLE]

Then, using Eq. (22) with $\xi_{0}=1/2$ and $\xi_{1}={D(\alpha-1)C_{0}K_{D}}/{4\left(K_{D}-\alpha C_{0}\right)^{2}}$, we obtain

[TABLE]

C.2 High-diffusion regime

Similar to the low $\tilde{D}$ case, we now assume the final solution for $\tilde{D}\gg 1$ includes a small first-order correction $\epsilon G_{1}$ to $G_{\textrm{opt}}$, but with $\epsilon\equiv 1/\tilde{D}$:

[TABLE]

We note here the optimal $G$ (high-$\tilde{D}$ limit) is constant. Again, with $G,G^{*}\gg K_{M}$, Eq. (23) gives, in the steady state:

[TABLE]

where in the first equality we have used $f_{u}=1-f_{b}$. Note that $f_{b}=C/(C+K_{\textrm{eff}})$ depends on the concentration profile, so we can expand it to first order in $g$ as $f_{b}(\tilde{G})\approx({1}/{2})+({g}/{8})\cos(\theta-\phi)$. Then, comparing the terms of order $\epsilon^{0}$ in Eq. (34):

[TABLE]

which can be simply solved via ansatz to give

[TABLE]

Finally, the correction is

[TABLE]

Following the same steps as Appendix C.1, we can get from (5):

[TABLE]

Then again, with $\xi_{0}=1/2$ and $\xi_{1}=\frac{1}{8}\left(1-\frac{\left(K_{D}-\alpha C_{0}\right)^{2}}{2\tilde{D}(\alpha-1)C_{0}K_{D}}\right)$, we arrive at the Fisher information:

[TABLE]

We can see this more easily as a small correction to the optimal Fisher information $Ng^{2}/32$.

Appendix D Simulation details

All simulations of the second-order PDEs like Eq. (13) are executed by Euler’s method of numerical integration using the forward-time centered-space method. We update (14) as follows:

[TABLE]

We update (40) for 40 points on the cell membrane ($\Delta\theta=\pi/20\approx 0.157$) until steady state and use Eq. (5) after fitting the bound probability to Eq. (19). Note that in doing so, we are assuming that the time scale associated with receptor-ligand binding is much shorter than the time scale associated with the allosteric protein’s activation/inactivation and diffusion.

For some parameter values, including for low values of $\tilde{D}$ and $C_{0}$, reaching steady state may require up to $\sim 3.3\times 10^{2}$ in unitless time. In these cases, we can both increase $\Delta t$ and the maximum number of time steps, the former of which can be afforded due to the low diffusion coefficient and the (approximate) stability condition $\frac{D\Delta t}{(\Delta x)^{2}}\lesssim\frac{1}{2}$. Unless we are sweeping over any of the parameters listed below, they remain at the values listed in Table 1. For the adaptation time calculations in Fig. 9, we run the simulation for a total time of $\tilde{T}=1000$ for the $1.2K_{D}/\alpha\to 1.1K_{D}/\alpha$ concentration jump, and $\tilde{T}=200$ for the other two. The concentration jump occurs at time $\tilde{t}_{0}=\tilde{T}/3$, allowing the system to reach steady state at both the initial and final concentrations. We then define the unitless adaptation time $\tilde{\tau}_{A}$:

[TABLE]

thus the physical $\tau_{A}$ can be computed from the unitless form as $\tau_{A}=\tilde{\tau}_{A}\times K_{G}/V_{\textrm{max}}$.