Universal properties of concentration sensing in large ligand-receptor networks

TL;DR

This paper demonstrates that large ligand-receptor networks can estimate multiple ligand concentrations from temporal binding sequences, revealing universal properties linked to Vandermonde matrices that can inform cellular biochemical decoding.

Contribution

It introduces a theoretical framework showing how cells can decode multiple ligand concentrations using receptor networks, highlighting universal spectral properties of the inverse covariance matrix.

Findings

Temporal sequences enable estimation of multiple ligands.

Universal spectral properties linked to Vandermonde matrices.

Potential implications for cellular biochemical decoding.

Abstract

Cells estimate concentrations of chemical ligands in their environment using a limited set of receptors. Recent work has shown that the temporal sequence of binding and unbinding events on just a single receptor can be used to estimate the concentrations of multiple ligands. Here, for a network of many ligands and many receptors, we show that such temporal sequences can be used to estimate the concentration of a few times as many ligand species as there are receptors. Crucially, we show that the spectrum of the inverse covariance matrix of these estimates has several universal properties, which we trace to properties of Vandermonde matrices. We argue that this can be used by cells in realistic biochemical decoding networks.

Figures

Equations

$$P(\{{\vec{\tau }}_i^{\mathrm{u}},{\vec{\tau }}_i^{\mathrm{b}}\}|\{c_{\alpha }\})$$

$$\times \sum _{\alpha =1}^{N_{\mathrm{L}}}{c}_{\alpha }{k}_{\alpha i}{r}_{\alpha i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\alpha i}}].$$

$$\begin{gathered} 0=\frac{\partial \log P}{\partial c_{\alpha }}\equiv (\log P{)}_{\alpha }^{\prime }=\frac{\partial \log \mathcal{P}}{\partial c_{\alpha }}+ \\ \sum _{i=1}^{N_{\mathrm{R}}}\sum _{m_i=1}^{n_i}\left(-{\tau }_{m_i}^{\mathrm{u}}{k}_{\alpha i}+\frac{k_{\alpha i}{r}_{\alpha i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\alpha i}}}{{\sum }_{\alpha =1}^{N_{\mathrm{L}}}{c}_{\alpha }{k}_{\alpha i}{r}_{\alpha i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\alpha i}}}\right), \end{gathered}$$

$$\log P^{\prime \prime }=\sum _{i=1}^{N_{\mathrm{R}}}\log P_i^{\prime \prime }+\log {\mathcal{P}}^{\prime \prime },$$

$${\left(\log P_i^{\prime \prime }\right)}_{\alpha \beta }=\sum _{m_i=1}^{n_i}\frac{(k_{\alpha i}{r}_{\alpha i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\alpha i}})(k_{\beta i}{r}_{\beta i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\beta i}})}{{\left({\sum }_{\alpha =1}^{N_{\mathrm{L}}}{c}_{\alpha }{k}_{\alpha i}{r}_{\alpha i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\alpha i}}\right)}^2}.$$

$$(J_i{)}_{m_i\alpha }=\frac{k_{\alpha i}{r}_{\alpha i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\alpha i}}}{{\sum }_{\gamma =1}^{N_l}{c}_{\gamma }{r}_{\gamma i}{e}^{-{\tau }_{m_i}^{\mathrm{b}}{r}_{\gamma i}}}.$$

$$(J_i{)}_{m_i\alpha }=\sum _{\nu =1}^{\infty }(J_i{)}_{\alpha }^{(\nu -1)}(d{\tau }_{m_i}{)}^{\nu -1},$$

$$\begin{gathered} ⟨({\tau }_i^{\mathrm{b}}{)}^n⟩=\sum _{\alpha }\frac{k_{\alpha i}{c}_{\alpha }}{{\sum }_{{\alpha }^{\prime }}{k}_{{\alpha }^{\prime }i}{c}_{{\alpha }^{\prime }}}\underset{0}{\overset{\infty }{\int }}{\tau }_{m_i}^n{r}_{\alpha i}{e}^{-r_{\alpha i}{\tau }_{m_i}}d{\tau }_{m_i} \\ =\sum _{\alpha }\frac{c_{\alpha }}{c_{\mathrm{tot}}}\frac{1}{r_{\alpha i}^n}\underset{0}{\overset{\infty }{\int }}{\tau }^n{e}^{-\tau }d\tau =\sum _{\alpha }\frac{c_{\alpha }}{c_{\mathrm{tot}}}\frac{\Gamma (n+1)}{r_{\alpha i}^n}, \end{gathered}$$

Full text

cellular information processing, random matrices, maximum likelihood, Vandermonde matrices

Introduction. Cellular environment contains many chemical ligands that are sensed by cell surface receptors. Typically, the number of ligand species is larger than that of the receptors. However, in traditional treatments of the problem, one only takes into account the (fluctuating) steady-state occupancy of receptors Berg and Purcell (1977); Bialek and Setayeshgar (2005); Hu et al. (2010); Kaizu et al. (2014), which allows estimation of just one quantity (e. g., one ligand species) per receptor. Recent work has focused instead on using sequences of binding and unbinding times for the estimation. The durations of unbound times carry information about the concentrations Endres and Wingreen (2009), while the durations of bound times identify the ligands. This allows to estimate the concentrations even in the presence of spurious background ligands Lalanne and Francois (2015); Mora (2015), or to get information about more than one concentration simultaneously from a single receptor Singh and Nemenman (2017); Carballo-Pacheco et al. (2019). All of this can be done using biologically plausible chemical reaction networks, such as variations of the kinetic proofreading model Hopfield (1974); Ninio (1975).

Here we turn to a previously not investigated regime, where both the number of ligands and the number of receptors are large, and ligands and receptors interact with a broad distribution of binding affinities, the so called Multiple Inputs – Multiple Outputs (MIMO) problem. We use maximum-likelihood (ML) techniques Endres and Wingreen (2009); Mora (2015); Singh and Nemenman (2017) to estimate concentrations of all ligands from outputs of all receptors, which provides estimates that are consistent with the true concentrations Wald (1949). We then focus on the co-variance of the estimates, obtained from the Hessian matrix. Our main finding is that the eigenvalue spectrum of this matrix exhibits universal behaviors, which we trace to properties of Vandermonde matrices. We argue that such MIMO problem is common in various biological systems, and it is also of relevance more broadly, beyond the ligand-receptor problem we study here.

Model. Consider a mixture of $N_{\rm L}$ ligands, with concentrations $c_{\rm\alpha},{\rm\alpha}\in[1,N_{\rm L}]$, that bind to $N_{\rm R}$ receptors with binding (unbinding) rates $k_{\alpha i}\;(r_{\alpha i})$. Notice that $N_{\rm L}$ can be larger than $N_{\rm R}$, so that more ligand concentrations are being measured than there are receptor types. Further, ligands and receptors cross-react, so that, in principle, $k_{\alpha i}>0$ and $r_{\alpha i}<\infty$ for all $\alpha$ and $i$. Thus we do not distinguish cognate and noncognate ligand-receptor pairs.

Suppose $n_{i}$ binding/unbinding events happen on the $i$’th receptor over the measurement time $T$. There is no way of knowing which particular ligand caused which binding, but binding durations are known, which provides some information about the ligand identity. In the regime of interest $n_{i}\gg 1$ and a single binding or unbinding event matters little. Thus we assume that all receptors are unbound at $t=0$ and bound at $t=T$. Then the sequence of unbound/bound durations of the $i$’th receptor is $\{\vec{\tau}_{i}^{\rm u},\vec{\tau}_{i}^{\rm b}\}=\{\tau^{\rm u}_{1_{i}},\tau^{\rm b}_{1_{i}},\tau^{\rm u}_{2_{i}},\tau^{\rm b}_{2_{i}},\dots,\tau^{\rm u}_{n_{i}},\tau^{\rm b}_{n_{i}}\}$. The likelihood of observing such a sequence is Singh and Nemenman (2017):

[TABLE]

The log likelihood can be maximized to get the maximum a posteriori (MAP) estimate of the concentration vector $\{c_{\alpha}\}$. Differentiating Eq. (1) w. r. t. $c_{\alpha}$ gives the following $N_{\rm L}$ coupled algebraic MAP equations:

[TABLE]

where $\mathcal{P}$ is the prior over the concentrations. The co-variance of the estimation can be obtained from the inverse of the Hessian matrix $(\log P^{\prime\prime})_{\alpha\beta}$ evaluated at the MAP solution. We get:

[TABLE]

where $\left(\log\mathcal{P}^{\prime\prime}\right)_{\alpha\beta}=\frac{\partial^{2}\log\mathcal{P}}{\partial c_{\alpha}\partial c_{\beta}}\delta_{\alpha\beta}$ and

[TABLE]

We assumed that the a priori concentration co-variances are zero. $\log P^{\prime\prime}_{i}$ contains terms that are contributed to the $\alpha\beta$ entry in the Hessian matrix by the receptor $i$. The sum in Eq. (4) is over the durations for the which this $i$’th receptor was bound, because only the bound durations are informative of the ligand identity.

Variances of the estimates scale as the inverse of the eigenvalues of $\log P^{\prime\prime}$ – high eigenvalues correspond to lower variance of various linear combination of the estimated concentrations. Thus we focus on the spectrum of eigenvalues of $\log P^{\prime\prime}$ in the rest of the paper.

Eigenvalue spectrum of the inverse covariance matrix. To illustrate the main properties of the eigenvalues of $\log P^{\prime\prime}$, we performed simulations of ligand-receptor networks of varying sizes. We explored different distributions of binding rates, concentrations, and unbinding rates, and these have little effect on the conclusions drawn below, unless noted otherwise. First, we set all (nominally unknown) ligand concentrations and binding rates to 1 for simplicity. This is the hardest inference problem – the least certainty about which specific ligand got bound to which specific receptor. Further, we choose to work with the log-normal distribution of unbinding rates (or, equivalently, normally distributed energy barrier between bound and unbound states). Figure 1 shows the eigenvalues of the matrix $\log P^{\prime\prime}$ vs the total simulation time $T$ for a network of 40 ligands and 4 receptors in panels (a) and (c), and a variable number of receptors in panel (b). In all cases, the highest singleton eigenvalue for each $T$ corresponds to the estimate of the total concentration, $c_{\rm tot}=\sum c_{\alpha}$. However, the rest of the eigenvalues come in subsets, whose size is equal to the number of receptors in the network, and these sets are separated almost equidistantly on a log-scale. Further, the mean of the eigenvalues and the top singleton eigenvalue are higher for systems with more receptors – which corresponds to a better inference coming from more independent samples of the concentrations, cf. Fig. 2. Finally, as the variance of the unbinding rates increases, the eigenvalues get lifted (the overall inference improves), and the splitting of subsets decreases, which suggests that (i) diversity of cross-reactivities among ligands and receptors improves the estimation, and (ii) the splitting is due to the degeneracy of nearly similar unbinding rates.

We now assume that the concentrations $c_{\alpha}$ are sampled from an unknown Gaussian distribution, $\mathcal{N}(\mu=1,\sigma=1)$. We can quantify the ability to infer the concentrations from the binding data by calculating the Fisher information – the expected value (over the concentration prior) of the Hessian of the log-likelihood matrix at the MAP solution. Specifically, we focus on the negative logarithm of the determinant of the Fisher information matrix, divided by the number of ligands, ${\mathcal{F}}\equiv-\langle\log P^{\prime\prime}_{\alpha\beta}|_{\rm MAP}\rangle/N_{\rm L}$. It represents the average logarithm of the variance on the ligand concentration following an observation. We investigate this quantity as a function of $N_{\rm R}$ and $N_{\rm L}$ in Figure 2. We observe that ${\mathcal{F}}\propto{N_{\rm L}}^{\alpha}/N^{\beta}_{\rm R}$, with $\alpha\approx 3.7$ and $\beta\approx 1$ for the observation duration $T=100$. The dependence on $N_{\rm R}$ represents the usual law of large numbers – as more measurements of the concentrations are made, the average variance of the estimates, ${\mathcal{F}}$, falls inversely proportionally with the number of measurements, $N_{\rm R}$. Correspondingly, $\beta$ does not change with $T$. The reason for scaling with $N_{\rm L}$ is unclear, but it indicates that decreasing the number of ligands has a dramatic effect on the ability to estimate them, approaching $\alpha\approx 4$ at asymptotically large $T$. Notice also that ${\mathcal{F}}$ starts deviating from 0 at $N_{\rm L}/N_{\rm R}\sim 3\dots 10$, indicating that the network can estimate nearly $3\dots 10$ as many ligand species as it has the receptors. Qualitatively similar results are obtained for other types distributions of the unbinding rates sup .

Origin of the eigenvalue spectrum. Two related arguments can explain our observations that (i) eigenvalues are split in groups of size $N_{\rm R}$, (ii) the groups are nearly equidistant from each other on the log space, (iii) higher diversity of unbinding rates decreases the gap between the groups. We start by writing $\log P^{\prime\prime}_{i}=J_{i}^{T}J_{i}$, where

[TABLE]

Assuming that the bound time durations on the receptor $i$, $\tau^{\rm b}_{m_{i}}$ are distributed narrowly around some mean value $\bar{\tau}^{\rm b}_{i}$ (the accuracy of this assumption depends on the tightness of the distribution of the unbinding rates), we can expand $(J_{i})_{m_{i}\alpha}$ around this mean value. This gives

[TABLE]

where $(J_{i})^{(\nu)}_{\alpha}$ represents the $\nu$’th derivative of $(J_{i})_{m_{i}\alpha}$, evaluated at $\bar{\tau}^{\rm b}_{i}$, and $d\tau_{m_{i}}=\tau_{m_{i}}^{\rm b}-\bar{\tau}^{\rm b}_{i}$,

Using the expansion in Eq. (6), we can write $J_{i}=V_{i}A_{i}$, where $\left(A_{i}\right)_{\nu\alpha}=(J_{i})^{(\nu-1)}_{\alpha}$ and $V_{i}$ is the Vandermonde matrix $(V_{i})_{m_{i}\nu}=\left(d\tau_{m_{i}}\right)^{\nu-1}$ Lam (1985). So, $\log P^{\prime\prime}_{i}=J_{i}^{T}J_{i}=(A_{i})^{T}(V_{i})^{T}V_{i}A_{i}$. The eigenvalues of the matrix $(V_{i})^{T}V_{i}$ scale as $d\tau^{2(\nu-1)}$ Waterfall et al. (2006). We performed simulations with matrices of this form, namely $(VR)^{T}(VR)$, where $R$ is a random matrix with elements chosen uniformly at random from $[0,1]$, and $d\tau_{i}$ in the Vandermonde matrices were chosen uniformly at random from $[-0.1,0.1]$. Other ranges and distributions of values produce qualitatively similar results. The resulting eigenvalues are shown in Fig. 3, together with eigenvalues of matrices $V$ and $R$, for comparison. These simulations suggest that the eigenvalues of $\log P^{\prime\prime}_{i}$ follow the same scaling as of $V^{T}V$.

The matrix formed by adding several matrices of the form $\log P^{\prime\prime}_{i}$ has an eigenvalue spectrum similar to that of $\log P^{\prime\prime}$ (compare Fig. 3 to Fig. 1). This is because the matrices being added ($\log P^{\prime\prime}_{i}$) have eigenvalues distributed roughly exponentially, but the corresponding eigenvectors are rotated randomly w. r. t. each other. In high-dimension, such random rotations result in the eigenvectors corresponding to the eigenvalue of the same rank being almost orthogonal to each other. This will introduce level splitting, similar to degenerate perturbation theory in quantum mechanics, so that if $N$ such matrices are added, eigenvalues will come in groups of $N$ sets.

We can illustrate the same result with a different, but related argument. We simulated matrices of the form $(M_{i})_{\alpha\beta}=\sum_{m=1}^{N}(x_{\alpha}x_{\beta})^{\tau_{m}}$, where $x_{\alpha}$’s and $\tau$’s are generated randomly, and $N\gg 1$ (our conclusions below hold for $N\gtrsim 5$). This form corresponds to the exponential terms in Eq. (4). The eigenvalues of these simulated matrices show an exponential scaling similar to that of $\log P^{\prime\prime}_{i}$, cf. Fig. 3. Further, if several such matrices are added together, $M=\sum_{i}M_{i}$, the resulting eigenvalue spectrum again looks similar to that of $\log P^{\prime\prime}$ because of the same level splitting argument for orthogonal eigenvectors corresponding to the same eigenvalues. This suggest that the eigenvalue spectrum we see for the ligand-receptor network is not overly specific to this system, but results from the particular structure (sum of exponentiated bilinear terms) of the Hessian matrices.

Decoding concentrations from receptor activities. Our receptor-ligand MIMO inference scheme only measures the duration of time for which each receptor is bound and unbound. Moments of these times can be used to infer concentrations of individual ligand species:

[TABLE]

where the term outside the integral in the first line is the probability of ligand $\alpha$ to bind receptor $i$ and $r_{\alpha i}e^{-r_{\alpha i}\tau_{m_{i}}}$ is the probability density for such bound interval durations. Here we assumed the binding rates to be the same for simplicity, which happens when they are diffusion-limited. Thus each moment is a linear combination of the ligand concentrations.

It is easy to design realistic biochemical networks that would solve such a systems of linear equations and infer the concentrations Mora (2015); Singh and Nemenman (2017); Carballo-Pacheco et al. (2019). For example, to estimate the first moment $\langle\tau_{i}\rangle$, a reporter molecule can be generated only when the receptor $i$ is bound. The mean reporter amount produced over time $T$ would be proportional to $T$, and, assuming many such molecules are produced over a typical bound interval, noise due to discreteness of the reporter would be negligible. Similarly, the estimate of the second moment, $\langle\tau_{i}^{2}\rangle$ can be obtained from a secondary reporter, which gets produced with the rate proportional to the instantaneous amount of the first reporter, and only while the receptor is bound. If the production rates are high, discreteness of this reporter will also be negligible. Then final network readouts can be activated / suppressed by the reporters to form their appropriate linear combinations representing $c_{\alpha}$Singh and Nemenman (2017).

Discussion. Here we studied a network of ligands and receptors with crosstalk, such that the number of ligands is larger than the number of receptors. Using the maximum likelihood solution based on time series of receptor binding and unbinding, we showed that estimation of ligands in this context, which could have been underdetermined, is, in fact, possible. We noticed that the maximum likelihood estimator results in a Hessian with a spectrum with a universal properties. Specifically, the eigenvalues come in subsets whose size is equal to the number of receptors in the network, and these subsets are almost exponentially distributed. This observation can be employed by biological systems to design simple and universal chemical kinetics schemes that would estimate concentrations of ligands from receptor activities irrespectively of the details of binding and unbinding rates in the network. Thus our theory makes specific predictions about the structure of molecular networks downstream of cell surface receptors. Whether biological systems follow these predictions remains to be seen, and networks around EGFR, IGF-1R, FC, and BMPR families of receptors Blinov et al. (2006); De Meyts and Whittaker (2002); Owen et al. (2009); Antebi et al. (2017), which respond to many related ligands by many related receptor, are good targets for the investigation. Finally, since we traced the universal spectrum of the Hessian to the properties of Vandermonde matrices that enter the inference problem, we expect that our findings will be applicable more generally, beyond ligand-receptor molecular networks, including problems in neuroscience and artificial sensory networks.

Acknowledgements. We thank Martin Tchernookov, Aditya Rangan, Thierry Mora and William S. Hlavacek for valuable discussions. VS was supported by University of Pennsylvania CNI fellowship. IN was supported in part by NIH and NSF Grants 1R01-EB022872, 5R01-NS099375, and PHY-1410978. IN acknowledges hospitality of the Kavli Institute for Theoretical Physics, supported in part by NSF and NIH Grants PHY-1748958, R25GM067110, and the Gordon and Betty Moore Foundation Grant 2919.01.