Introduction

A fundamental way of learning about a material is by observing how it responds to external stimuli. The functional dependence of a response on a stimulus is known as a constitutive relation. The most basic example of such a relation is Hooke’s law F = kX for the deformation response X of a linear elastic solid to a force stimulus F, where k is a material constant that is a characteristic property of the solid1,2. This linearity is a generic feature of material response for small enough stimuli, as it requires only that the constitutive relation be analytic and non-vanishing to first order. Linear constitutive relations have proven useful for characterizing a broad range of physical systems, including dielectric materials3, diffusion4, friction5, geomaterials6, Newtonian fluids7, piezoelectric materials8, thermoelectric materials9, and even abstract entities such as financial markets10,11.

Material constants of linear constitutive relations are typically inferred by comparing the known value of an applied stimulus to the measured response produced by the stimulus. For the case of a homogeneous elastic solid, the material constant is simply given by k = F/X. In reality, however, all materials are heterogeneous on a small enough scales12,13,14,15. This heterogeneity serves as a source of measurement noise that becomes significant for systems that operate at the microscale, such as miniature electronic devices16,17,18,19, medical microrobots20,21,22,23, and biological sensors24,25,26,27,28,29,30.

Previous studies of sensing in random media have focused on remote sensing or communication via traveling waves31,32,33,34,35,36. The inference of material properties at small scales has been studied in microrheology37,38,39 and for chemical sensing40,41,42,43. In these contexts, the measurement noise due to thermal fluctuations has been characterized using fluctuation–dissipation theorems41,44. However, thermal fluctuations are fundamentally different from structural heterogeneities: the former arise uniformly in space and vary in time, and the latter vice versa. Although methods are available to probe heterogeneous materials on small scales, it is not known how precisely this process can be done13,45,46,47,48. What are the limits to sensing the properties of heterogeneous materials, and how can a physical device be designed to achieve these limits?

To quantify the limits of sensing constitutive relations, we investigate a simple model of a localized sensor that probes a heterogeneous medium to estimate a global material constant. Specifically, we consider a continuous medium with a material constant given by a uniform average value λ0 plus a spatially varying fluctuation δλ(r) with short-ranged correlations. We treat the sensor as a spherical device that can probe λ0 by applying an external stimulus field and measuring the resulting response field in equilibrium.

In what follows, we show that this inference process admits an optimal (minimum-variance unbiased) measurement protocol. Surprisingly, the optimal protocol depends qualitatively on both the spatial resolution of the sensor and also on whether it can perform multiple probes. For a single probe, the optimal protocol is remarkably complex, because the modes applied to the medium can interfere with each other in a geometrically frustrated manner, akin to the spins of a spin-glass. In contrast, the optimal protocol for multiple probes is comparatively simpler, because it avoids unnecessary interference effects. We exploit this simplicity to determine the total amount of information that the sensor can extract by probing a given region of the medium. Physically, optimal performance is achieved by decoding the results of a sequence of probes that penetrate into the surrounding medium to varying extents. This strategy can allow the sensor to effectively average the material constant over a volume that grows with the spatial resolution of the sensor. Finally, we use our theoretical framework to bound the precision of mechanosensing in a biopolymer network, a sensory process that regulates cellular behavior in decisive ways28,49,50,51 and is used for medical diagnostics and material fabrication52,53,54.

Results

Probing a Winkler foundation

To gain insight into sensing material properties in physical space, we explore a minimal model that consists of a spherical sensor embedded in a heterogeneous medium (Fig. 1). In this section, we start by taking the medium to be the simplest heterogeneous material: a disordered Winkler foundation55. This medium corresponds to an array of decoupled springs in the continuum limit. The internal energy of the Winkler foundation is given by:

$$E=\frac{1}{2}\int\lambda ({\boldsymbol{r}})u{({\boldsymbol{r}})}^{2}{\mathrm{d}}{\boldsymbol{r}},$$
(1)

where λ(r) is a spatially varying material constant and u(r) is the response field at position r. We assume λ(r) = λ0 + δλ(r), where λ0 is a fixed, uniform field and δλ(r) λ0 is a Gaussian random field with zero mean and spatial correlations given by:

$$\langle \delta \lambda ({\boldsymbol{r}})\delta \lambda ({\boldsymbol{r}}^{\prime} )\rangle =\frac{{\Delta }_{\lambda }}{{(2\pi )}^{D/2}}{e}^{-{({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} )}^{2}/{\xi }^{2}},$$
(2)

where \({\Delta }_{\lambda }\ll {\lambda }_{0}^{2}\) is the local variability of λ(r), D is the spatial dimension, and ξ is the correlation length of the fluctuations in λ(r). For simplicity, we assume ξ is small enough that these correlations can be approximated by:

$$\langle \delta \lambda ({\boldsymbol{r}})\delta \lambda ({\boldsymbol{r}}^{\prime} )\rangle ={\Delta }_{\lambda }{\xi }^{D}\delta ({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} ).$$
(3)
Fig. 1: Sensing in a heterogeneous medium.
figure 1

Schematic illustration of sensing model, showing an idealized, spherical sensor of radius a (green) embedded inside a medium with a spatially varying material constant field λ(r) (background). The sensor can learn about λ(r) by applying an arbitrary stimulus and recording an arbitrary weighted response within its volume.

The quenched disorder δλ(r) in the material constant limits the precision with which a physical sensor can infer λ0. To determine these limits, we consider an idealized sensor that probes λ0 by first applying a stimulus field f(r). This field perturbs the energy of the system as follows:

$$\delta E=-\int\,f({\boldsymbol{r}})u({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}.$$
(4)

After applying this stimulus, the sensor measures the response of the medium in equilibrium. In particular, we assume that the sensor records an integrated response m:

$$m=\int\ w({\boldsymbol{r}})u({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}},$$
(5)

where w(r) is a weight field. Taken together, the probe fields f(r) and w(r) define the measurement protocol of the sensor. For any physical sensor, these fields must be localized in space. We impose this locality by constraining the probe fields to obey f(r) = 0 and w(r) = 0 for r > a, where r is the radial coordinate and a is the radius of the sensor.

Finally, upon recording the integrated response m, the sensor produces an estimate for λ0. In what follows, we will determine the optimal estimator \({\hat{\lambda }}_{0}\) perturbatively to leading order in δλ(r). In this approximation, the integrated response is:

$$m=\int\ \left(\frac{1}{{\lambda }_{0}}-\frac{\delta \lambda ({\boldsymbol{r}})}{{\lambda }_{0}^{2}}\right)\psi ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}},$$
(6)

where we have defined the probe intensity ψ(r) ≡ f(r)w(r), a function of the probe fields that captures how strongly the probe senses a given location of space. For a fixed choice of ψ(r), along with prior knowledge of the model parameters other than λ0, the optimal estimator of λ0 based on the outcome of m is (Supplemental Material, Supplementary Note 1):

$${\hat{\lambda }}_{0}=\frac{s}{m},$$
(7)

where s is a normalizing constant chosen such that the estimator \({\hat{\lambda }}_{0}\) yields an unbiased estimate of λ0:

$$s=\int\ \psi ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}.$$
(8)

Equation (7) is a mesoscopic generalization of Hooke’s law k = F/X. By computing the estimate \({\hat{\lambda }}_{0}\), the sensor obtains a weighted spatial average of λ(r):

$${\hat{\lambda }}_{0}=\frac{\int\ \psi ({\boldsymbol{r}})\lambda ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}}{\int\ \psi ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}},$$
(9)

to leading order in δλ(r). This estimator is optimal in that it has a lower standard deviation \(\delta {\lambda }_{0}\equiv \sqrt{\langle {({\hat{\lambda }}_{0}-{\lambda }_{0})}^{2}\rangle }\) than any other unbiased estimator for a fixed choice of ψ(r). Therefore, the optimal measurement protocol can be determined by minimizing \(\delta {\lambda }_{0}^{2}\) with respect to the probe intensity ψ(r). Inserting Eq. (3) into the definition of the variance yields:

$$\delta {\lambda }_{0}^{2}={\Delta }_{\lambda }{\xi }^{D}\frac{\int\ \psi {\left({\boldsymbol{r}}\right)}^{2}{\mathrm{d}}{\boldsymbol{r}}}{{\left(\int\psi ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}\right)}^{2}}.$$
(10)

This variance is invariant with respect to an overall rescaling of ψ(r). To eliminate this redundancy, we constrain ∫ ψ(r)dr to be a fixed constant. Furthermore, we must enforce ψ(r) = 0 in the exterior of the sensor (r > a) to satisfy the constraints imposed by the finite size of the sensor. Thus, the minimum of \(\delta {\lambda }_{0}^{2}\) is determined by the configuration of ψ(r) that extremizes the following action S:

$$S={\int}_{{{\mathcal{R}}}_{{\rm{int}}}}\left(\frac{1}{2}\psi {\left({\boldsymbol{r}}\right)}^{2}-\gamma \psi ({\boldsymbol{r}})\right){\mathrm{d}}{\boldsymbol{r}},$$
(11)

where the integral is taken over the interior \({{\mathcal{R}}}_{{\rm{int}}}\) of the sensor (r < a) and γ is a Lagrange multiplier that fixes ∫ψ(r)dr. This action is extremized by any measurement protocol with a probe intensity ψ(r) that is uniform over \({{\mathcal{R}}}_{{\rm{int}}}\). The optimal measurement protocol is, therefore:

$$\psi ({\boldsymbol{r}})=\left\{\begin{array}{ll}\gamma ,&r\, <\, a.\\ 0,&r \,> \, a.\end{array}\right.$$
(12)

Inserting Eq. (12) into Eq. (10) yields:

$$\delta {\lambda }_{0}^{2}={\Delta }_{\lambda }{\xi }^{D}{V}^{-1},$$
(13)

where V is the volume of the sensor. Thus, the fractional uncertainty of the estimator \({\hat{\lambda }}_{0}\), defined as the standard deviation δλ0 divided by the mean λ0, scales as:

$$\frac{\delta {\lambda }_{0}}{{\lambda }_{0}} \sim {\left(\frac{{\Delta }_{\lambda }}{{\lambda }_{0}^{2}}\right)}^{1/2}{\left(\frac{\xi }{a}\right)}^{D/2},$$
(14)

which can be interpreted as the familiar \(1/\sqrt{N}\) scaling of measurement uncertainty for N independent samples. In this analogy, the sample size N ~ (a/ξ)D corresponds to the number of effectively independent subvolumes probed by the sensor.

Probing an elastic sheet

For the Winkler foundation, our model sensor could not induce a response beyond its volume. In contrast, many other types of elastic media are coupled in space and thereby respond to stimuli nonlocally. To understand how such nonlocality affects a sensor’s ability to infer material properties, we now turn to conventional, linear elasticity. For simplicity, we will first focus on an isotropic, two-dimensional elastic sheet characterized by a single material constant, and in section “The precision of biomechanical sensing”, we will generalize our theoretical framework to a three-dimensional elastic medium characterized by a material constant tensor.

For the elastic sheet, we consider the deformation response u(r) to force stimuli f(r) oriented perpendicular to the plane of the sheet. Thus, the sheet’s internal energy depends on the gradient u(r) of the response field as follows:

$$E=\frac{1}{2}\int\ \lambda ({\boldsymbol{r}})\nabla u({\boldsymbol{r}})\cdot \nabla u({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}.$$
(15)

Here, as in the previous section, we take λ(r) to be a Gaussian random field with mean λ0, variance \({\Delta }_{\lambda }\ll {\lambda }_{0}^{2}\), and spatial correlations over a scale ξ. As before, we take the sensor to interact with the medium within a radius a by first applying a stimulus field f(r) as in Eq. (4), and then measuring an integrated response m as in Eq. (5).

To leading order in δλ(r), the sensor can again compute \({\hat{\lambda }}_{0}=s/m\) to obtain a spatial average of λ(r) weighted by a probe intensity ψ(r), as in Eq. (9) (Supplemental Material, Supplementary Note 2). However, for the elastic sheet, ψ(r) is now:

$$\psi ({\boldsymbol{r}})={\left(\nabla \cdot \right)}^{-1}f({\boldsymbol{r}})\cdot {\left(\nabla \cdot \right)}^{-1}w({\boldsymbol{r}}),$$
(16)

where ()−1 is the inverse divergence operator (Supplementary Note 2). This probe intensity is a nonlocal function of the probe fields and thereby allows the sensor to probe distant regions beyond its boundary.

Intuitively, probing a greater extent of the medium should yield a more accurate estimate of λ0. To that end, the greatest possible extent of a probe is achieved by probe potentials with a ~1/r radial dependence in the far-field limit. For the elastic sheet, this decay profile is not produced by monopoles (which yield pathological, non-decaying potentials), but rather by dipoles. The simplest possible measurement protocol with dipole probe fields is described by:

$$f({\boldsymbol{r}}) \sim \delta (r-a)\cos (\theta ),$$
(17)
$$w({\boldsymbol{r}}) \sim \delta (r-a)\cos (\theta ).$$
(18)

These probe fields cast a probe intensity ψ(r) that is uniform in the interior of the sensor and isotropically decaying in the exterior:

$$\psi ({\boldsymbol{r}})=\left\{\begin{array}{ll}\gamma ,&r \,< \,a.\\ \gamma {\left(\frac{a}{r}\right)}^{4},&r\, > \,a.\end{array}\right.$$
(19)

Inserting Eq. (19) into Supplementary Equation 31 yields the following variance:

$$\delta {\lambda }_{0}^{2}=\frac{1}{3}{\Delta }_{\lambda }{\xi }^{D}{V}^{-1},$$
(20)

for D = 2. As expected from dimensional analysis, this expression has the same dependence on the model parameters as for the Winkler foundation (cf. Eq. (13)). Importantly, however, its prefactor is smaller. Thus, our example illustrates how a sensor can harness a long-ranged response function to perform at a higher precision by effectively averaging λ(r) over a larger region of space.

Probe-field interference limits the channel capacity of sensing

Given that a probe of the elastic sheet can access nonlocal information, what limits its precision? To answer this question, we start by considering the simpler case of a sensor that can only apply probe fields on its boundary. For such boundary probes, the most general probe fields are of the form:

$$f({\boldsymbol{r}}) \sim \delta (r-a)\sum _{k}{B}_{k}^{(f)}{e}^{ik\theta },$$
(21)
$$w({\boldsymbol{r}}) \sim \delta (r-a)\sum _{k}{B}_{k}^{(w)}{e}^{-ik\theta },$$
(22)

where \({B}_{k}^{(f)}\) and \({B}_{k}^{(w)}\) are complex coefficients that satisfy \({B}_{-k}^{(f)}={B}_{k}^{(f)* }\) and \({B}_{-k}^{(w)}={B}_{k}^{(w)* }\) to ensure that the probe fields are real, and we assume k > 0 to avoid pathological, non-decaying interactions caused by monopoles. This measurement protocol casts the following probe intensity:

$${\psi }_{\pm }({\boldsymbol{r}})=\sum _{k,l}{B}_{kl}(kl+| k| | l| ){\left(\frac{r}{a}\right)}^{\pm | k| \pm | l| -2}{e}^{i(k-l)\theta },$$
(23)

where ψ+(r) and ψ(r) are probe intensities that correspond to the interior (r < a) and the exterior (r > a) of the sensor, respectively, and \({B}_{kl} \sim {B}_{k}^{(f)}{B}_{l}^{(w)}\). Inserting ψ±(r) into the definition of the variance and performing the spatial integrals yields:

$$\delta {\lambda }_{0}^{2}={\Delta }_{\lambda }{\xi }^{D}\sum _{k,l,m,n}{B}_{kl}{B}_{mn}{{\mathbb{T}}}_{klmn},$$
(24)

where \({{\mathbb{T}}}_{klmn}\) is a highly structured, fourth-order tensor:

$${{\mathbb{T}}}_{klmn}=4\pi {a}^{2}\frac{{\delta }_{k-l+m-n,0}{x}_{klmn}{y}_{klmn}}{({x}_{klmn}+2)({x}_{klmn}-2)}.$$
(25)

Here, δi,j is the Kronecker delta function, \({x}_{klmn}=\left|k\right|+\left|l\right|+\left|m\right|+\left|n\right|\), and \({y}_{klmn}=(kl+\left|kl\right|)(mn+\left|mn\right|)\). In Eq. (24), we have normalized ψ±(r) such that ∫ψ(r)±dr = 1, which implies that Bkl must obey:

$$\sum _{k}4\pi {a}^{2}| k| {B}_{kk}=1.$$
(26)

To gain insight into the optimal measurement protocols for boundary probes, we used the Nelder–Mead algorithm to numerically minimize δλ0/λ0 over \({B}_{k}^{(f)}\) and \({B}_{k}^{(w)}\) (see “Methods” section). To that end, we imposed a cutoff on the system by truncating the sums in Eqs. (24) and (26) at a maximum absolute mode number \({k}_{\max }\). Physically, this parameter corresponds to the spatial resolution of the sensor, which we define as:

$$d=\left(\frac{2\pi }{{k}_{\max }}\right)a.$$
(27)

For all choices of \({k}_{\max }\), we studied, this algorithm converged to basins of minima dominated by the dipole modes (k = 1), which makes intuitive sense given that these modes probe the largest extent of the medium. Interestingly, however, as we increased \({k}_{\max }\), we found that at certain special values, the optimal probe fields shifted and picked up additional higher-order modes, resulting in a smaller minimum fractional uncertainty \(\delta {\lambda }_{0,\min }/{\lambda }_{0}\) (Fig. 2a).

Fig. 2: Probe-field interference can limit the information that a sensor can glean from a single probe.
figure 2

a Fractional uncertainty \(\delta {\lambda }_{0,\min }/{\lambda }_{0}\) in units of η = ξDΔλV−1 for numerically optimal boundary probes versus maximum absolute mode number \({k}_{\max }\) included in the multipole expansions of the probe fields (see section “Probe-field interference limits the channel capacity of sensing”). Red point corresponds to \({k}_{\max }=1\), green point to \({k}_{\max }=2\), and blue points to \({k}_{\max }> \,2\). Dashed gray line indicates the optimum \(\delta {\lambda }_{0,{\rm{low}}}/{\lambda }_{0}\approx \eta /\sqrt{\pi }\) attained in the limit \({k}_{\max }\to \infty\) for the convex relaxation of δλ0/λ0 described in section “Probe-field interference limits the channel capacity of sensing”. Inset: \(\delta {\lambda }_{0,\min }/{\lambda }_{0}-\eta /\sqrt{\pi }\) in units of η versus \({k}_{\max }\) on a logarithmic scale. b Absolute values of the weight field coefficients \(| {B}_{k}^{(w)}|\) versus the absolute values of the stimulus field coefficients \(| {B}_{k}^{(f)}|\) for an example measurement protocol obtained via numerical optimization for \({k}_{\max }=16\), showing dipole modes (red), quadrupole modes (green), and higher-order modes (blue). Dashed gray line shows \(| {B}_{k}^{(w)}| =| {B}_{k}^{(f)}|\). Inset: phases of \({B}_{k}^{(f)}\) and \({B}_{k}^{(w)}\) for the same measurement protocol as in the main panel. Lines connect the coefficients that correspond to the same value of k. Colors same as in the main panel. c Probe intensity ψ(r) for the same measurement protocol as in b versus spatial coordinate r. Left inset: ψ(r) at the boundary of the sensor (r = a) versus angular coordinate θ. Right inset: larger view of the region indicated by the black rectangle in the main panel, showing small wrinkles in ψ(r).

The higher-order modes contribute with smaller amplitudes and nontrivial relative phase shifts (Fig. 2b). These complex configurations arise because different terms in Eq. (24) can provide conflicting contributions to \(\delta {\lambda }_{0}^{2}\) depending on the relative phases of the modes. This geometrical frustration greatly suppresses modes beyond the dipole-dipole and quadrupole-quadrupole pairs, which for \(2\le {k}_{\max }\,\le \,12\) appear together with amplitudes and phase relations that maximize the extent of ψ(r) while preserving its isotropy. Including three or more mode pairs must break isotropy, analogous to how three or more antiferromagnetic spins cannot simultaneously minimize their interaction energies (Supplementary Note 3). Nevertheless, for \({k}_{\max }\,> \,12\), the optimal measurement protocols contain additional higher-order modes that cause small wrinkles in ψ(r) (Fig. 2c). Although these wrinkles break the isotropy of ψ(r), they also smoothen out its profile in the radial direction, which results in a greater overall uniformity throughout space and thus a higher precision.

To better understand the asymptotic behavior of \(\delta {\lambda }_{0,\min }/{\lambda }_{0}\) for large \({k}_{\max }\), we imagine relaxing the constraints on \(\delta {\lambda }_{0}^{2}\) by allowing Bkl to be an arbitrary matrix satisfying \({B}_{-k,-l}={B}_{kl}^{* }\). This relaxation expands the space of possible ψ(r) to include all real configurations that can be generated by Eq. (23), some of which cannot be cast by a physical probe. Importantly, this relaxation is a convex function of Bkl, and thus has a unique minimum δλ0,low/λ0 that provides a theoretical lower bound on δλ0/λ0. Specifically, in the limit \({k}_{\max }\to \infty\), we find that \(\delta {\lambda }_{0,{\rm{low}}}/{\lambda }_{0}\approx {\xi }^{D}{\Delta }_{\lambda }{V}^{-1}/\sqrt{\pi }\), which provides a close lower bound on the values obtained via numerical minimization (Fig. 2a and Supplementary Note 4).

A simple argument based on symmetry reveals that this lower bound must be a strict inequality for \({k}_{\max }> 2\). This argument follows from observing that for all \({k}_{\max }\), the unique optimal configuration of ψ(r) for the relaxation is isotropic, in contrast to the configurations we found by minimizing Eq. (24) for \({k}_{\max }\,> \,12\) (Supplementary Note 4). This broken isotropy must persist for all higher values of \({k}_{\max }\), and therefore a boundary probe can never cast a configuration of ψ(r) that performs as well as the optimal ψ(r) for the convex relaxation of δλ0/λ0. This example illustrates how interferences between the probe fields limit the information that can be gleaned from a single probe, i.e. the channel capacity of sensing. In the following section, we will show how a sensor can overcome this limit by performing multiple probes, and then we will generalize our results to a sensor that can apply arbitrary probe fields in its volume.

Sensory multiplexing can significantly improve the precision of sensing

The interferences in the previous section occur because all of the modes applied by a probe interrogate the medium simultaneously. In principle, however, each mode couples to a different spatial extent of the medium and therefore should carry independent information about λ0. Such information could potentially be accessed by performing separate measurements with distinct spectra.

To test this notion, we determine the optimal estimator for a sensor that can perform multiple probes with varying measurement protocols. We label each probe by an integer k and constrain their probe fields to be zero for r > a. In this case, the minimum-variance unbiased estimator of λ0 is again given by a weighted spatial average of λ(r):

$${\hat{\lambda }}_{0}=\frac{\int\ \Psi ({\boldsymbol{r}})\lambda ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}}{\int\ \Psi ({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}}.$$
(28)

Here, Ψ(r) is an effective probe intensity created by the optimally weighted sum of the probe intensities ψk(r) for the individual probes:

$$\Psi ({\boldsymbol{r}})=\sum _{k}{p}_{k}\frac{{\psi }_{k}({\boldsymbol{r}})}{\int\ {\psi }_{k}({\boldsymbol{r}}){\mathrm{d}}{\boldsymbol{r}}},$$
(29)

where \({p}_{k}={\sum }_{l}{C}_{kl}^{-1}\) with \({C}_{kl}\equiv \langle ({\hat{\lambda }}_{0,k}-{\lambda }_{0})({\hat{\lambda }}_{0,l}-{\lambda }_{0})\rangle\) defined as the covariance matrix of the estimators \({\hat{\lambda }}_{0,k}\) for the individual probes (Supplementary Note 5). The variance of \({\hat{\lambda }}_{0}\) is:

$$\delta {\lambda }_{0}^{2}={\left(\sum _{k,l}{C}_{kl}^{-1}\right)}^{-1}.$$
(30)

For simplicity, we start by considering a sensor that applies a sequence of probe fields:

$${f}_{k}({\boldsymbol{r}}) \sim \delta (r-a)\cos (k\theta ),$$
(31)
$${w}_{k}({\boldsymbol{r}}) \sim \delta (r-a)\cos (k\theta ).$$
(32)

from an initial mode number k = 1 up to a maximum mode number \(k={k}_{\max }\), which corresponds to the spatial resolution d of the sensor defined by Eq. (27). By varying k, the sensor modulates the range of ψk(r) in the exterior at the cost of simultaneously modulating ψk(r) in the interior:

$${\psi }_{k}({\boldsymbol{r}}) \sim \left\{\begin{array}{ll}{\left(\frac{r}{a}\right)}^{2k-2},&r\,< \,a.\\ {\left(\frac{r}{a}\right)}^{-2k-2},&r\,> \, a.\end{array}\right.$$
(33)

Interestingly, this collection of boundary probes does not achieve a significant improvement over a single optimal boundary probe. Instead, as \({k}_{\max }\) is increased, the fractional uncertainty approaches \(\delta {\lambda }_{0}/{\lambda }_{0}\approx {\xi }^{D}{\Delta }_{\lambda }{V}^{-1}/\sqrt{\pi }\), as we found for the convex relaxation in section “Probe-field interference limits the channel capacity of sensing”. This agreement is not a mere coincidence: for boundary probes, the possible configurations of Ψ(r) are mathematically equivalent to the possible configurations of ψ(r) for the convex relaxation of a single probe (Supplementary Note 6). However, unlike the convex relaxation, the collection of boundary probes reveals an additional physical effect that can limit the precision of a sensor. That is, for multiple probes, the overlapping configurations of ψk(r) in the interior correlate the probes and thereby suppress the amount of information that can be extracted from the exterior. These correlations are reflected in the structure of the covariance matrix:

$${C}_{kl}=\frac{1}{4}{\Delta }_{\lambda }{\xi }^{2}{V}^{-1}\left(\frac{kl}{k+l-1}+\frac{kl}{k+l+1}\right).$$
(34)

In this expression, the first and second fractions are contributed by overlaps in the interior and exterior, respectively. To compensate for the superfluous contributions from the interior, the sensor must employ probe fields that are nonzero within its volume. One way to perform this compensation is by pairing each probe k with a companion probe described by (Supplementary Note 6):

$${\tilde{\psi }}_{k}({\boldsymbol{r}}) \sim \left\{\begin{array}{ll}{\left(\frac{r}{a}\right)}^{2k-2},&r \,< \,a.\\ 0,&r \,> \,a.\end{array}\right.$$
(35)

Pairing these companion probes with the original probes using Eq. (29) with appropriate values of pk yields effective probe intensities Ψk(r) that are zero in the interior (Supplementary Note 6):

$${\Psi }_{k}({\boldsymbol{r}}) \sim \left\{\begin{array}{ll}0,&r \,< \,a.\\ {\left(\frac{r}{a}\right)}^{-2k-2},&r \,> \,a.\end{array}\right.$$
(36)

Finally, the sensor may include an additional unpaired companion probe \({\tilde{\psi }}_{k}({\boldsymbol{r}})\) with k = 1 to uniformly sample λ(r) in its interior. With these adjustments, the resulting all-inclusive effective probe intensity Ψ(r) exhaustively decodes the information that can be extracted by probing the region r < a using any combinations of probe fields (Supplementary Note 7). In this case, the covariances among the paired probes Ψk(r) and the unpaired probe \({\tilde{\psi }}_{1}({\boldsymbol{r}})\) are given by:

$${C}_{kl}={\Delta }_{\lambda }{\xi }^{2}{V}^{-1}\left({\delta }_{k,0}{\delta }_{l,0}+\frac{kl}{k+l+1}\right),$$
(37)

for kl ≥ 0, where kl = 0 correspond to the unpaired probe. Inserting the inverse of this matrix into Eq. (30) yields:

$$\delta {\lambda }_{0}^{2}={\Delta }_{\lambda }{\xi }^{2}{V}^{-1}{\left(\frac{1}{{k}_{\max }+1}\right)}^{2}.$$
(38)

This variance decreases with \({k}_{\max }\) because each additional probe increases the uniformity of Ψ(r) over space (Fig. 3). In the limit of fine resolution \({k}_{\max }\gg 1\) (d a), the fractional uncertainty of the sensor’s estimate scales as:

$$\frac{\delta {\lambda }_{0}}{{\lambda }_{0}} \sim {\left(\frac{{\Delta }_{\lambda }}{{\lambda }_{0}^{2}}\right)}^{1/2}{\left(\frac{d}{a}\right)}^{D/2}{\left(\frac{\xi }{a}\right)}^{D/2},$$
(39)

for D = 2. Thus, simultaneously varying both probe fields throughout the volume of the sensor can allow a significant amount of additional information to be transmitted across the sensory channel. We refer to the strategy of performing multiple measurements with varying probe fields as “sensory multiplexing.”

Fig. 3: Sensory multiplexing can greatly improve the precision of sensing.
figure 3

a Colored points show the smallest attainable fractional uncertainty δλ0/λ0 in units of η = ξDΔλV−1 for a sensor that can perform sensory multiplexing up to a maximum absolute mode number \({k}_{\max }\) (see section “Sensory multiplexing can significantly improve the precision of sensing”). Gray circles show a lower bound δλ0,low/λ0 on the fractional uncertainty for each value of \({k}_{\max }\) for a single volume probe, obtained via numerical minimization (Supplementary Note 10). b All-inclusive effective probe intensities Ψ(r) for sensory multiplexing versus radial coordinate r in units of the sensor radius a for the same values of \({k}_{\max }\) as in a (correspondence indicated by matching colors). Inset shows Ψ(r) on a lin-log scale.

Sensory multiplexing can be generalized to D = 3 by taking the probe fields to be pairs of spherical harmonics. In this case, δλ0/λ0 still obeys the asymptotic scaling in Eq. (39) (Supplementary Note 8). Moreover, this scaling is robust to the omission of a finite number of modes (Supplementary Note 9). Taken together, our results reveal that for d a, sensory multiplexing can improve the fractional uncertainty of a sensor by a factor proportional to the number (a/d)D of distinct subvolumes that it can resolve simultaneously for D = 2 and D = 3.

Notably, this level of precision can never be attained by a single probe, even if the sensor is permitted to apply an arbitrary pair of probe fields within its volume. This limitation occurs due to probe-field interference, as before in section “Probe-field interference limits the channel capacity of sensing”. That is, including more than three pairs of boundary modes breaks the isotropy of ψ(r), and a sensor can always improve upon an anisotropic ψ(r) by performing multiple rotated copies of the probe and combining the results using Eq. (28). Moreover, numerical optimization suggests that Eq. (38) does not provide a close bound on the precision of a single volume probe, even if the sensor is allowed to separately optimize ψ(r) in the interior and the exterior (Fig. 3a and Supplementary Note 10). In sum, we conclude that the optimal sensing strategy involves combining the results of multiple different measurements that probe the medium with different multipole symmetries.

The precision of biomechanical sensing

In this section, we extend our modeling framework to a scenario in which structural heterogeneity is known to play a significant role: cellular mechanosensing. Certain types of eukaryotic cells actively probe and respond to the stiffness of their surroundings, which has been shown to guide their behavior in decisive ways24,56,57,58,59,60. The importance of such mechanosensing invites the question of how precisely cells exploit the mechanical information available to them. In what follows, we present numerical evidence suggesting that some cells make optimal use of this information.

In connective tissue, a cell’s mechanical environment primarily consists of a disordered biopolymer network that serves as a scaffold on which the cell lives and moves28,29,61,62. To quantify what a cell can learn by interacting with such a network, we generalize our sensing model to a three-dimensional, isotropic elastic medium characterized by a shear modulus μ and a Poisson’s ratio σ (Supplementary Note 11). For simplicity, we take σ to be a fixed, uniform field and μ to be the sum of a fixed, uniform field μ0 and a spatially-varying random field δμ(r) with short-ranged correlations as in Eq. (3).

We determined the parameters in our model for a reconstituted collagen network, an in vitro system that closely resembles in vivo cellular environments28,29,49,50. For a collagen network prepared from a c ~ 0.2 μg/mL solution of collagen type-I monomers, previous studies suggest μ0 ~ 0.3 Pa, σ ~ 0.4, Δμ ~ 0.1 Pa2, and ξ ~ 5 μm (Supplementary Note 11). For these values, the ratio \({\Delta }_{\mu }^{1/2}/{\mu }_{0} \sim 1\) lies outside the strict regime of validity of our perturbative approach; nevertheless, we expect our results to qualitatively describe how δμ0/μ0 depends on the model parameters to the right order of magnitude.

Eukaryotic cells can sense stiffness by attaching to biopolymer networks via transmembrane protein complexes called focal adhesions28,56,63,64. For D = 3, the simplest cellular probe consists of isotropic dipolar shells of radius a (Supplementary Note 11). Taking a = 10 μm leads to δμ0/μ0 ~ 0.15, which supports the notion that cells could use mechanical information to reliably distinguish between different connective tissue environments, including brain (μ0 ~ 1 kPa), muscle (μ0 ~ 10 kPa), and bone (μ0 ~ 100 kPa)28,57,58,59,60. Such mechanosensing could be tested in experiment by using micropatterned materials to explore the effect of substrate heterogeneity on intracellular signaling dynamics. Interestingly, previous work in which cells were seeded on two-dimensional patterned substrates found that increases in either the average stiffness or in the spatial uniformity of substrate heterogeneities both resulted in increased intracellular signal transduction27. This agreement is consistent with Eq. (39), which suggests that increasing μ0 and decreasing ξ should both have the same sign of influence on cellular precision.

A cell could further reduce δμ0/μ0 via sensory multiplexing. Indeed, cells have been known to modulate the forces they exert in order to vary the modes they apply, consistent with our predictions for behavior during sensory multiplexing65. A cell’s spatial resolution is limited by the maximum number of focal adhesions that it can simultaneously apply to the network. Cells have been observed to display more than  ~100 focal adhesions66, which could allow a 10 μm cell to probe the network on scales smaller than ξ. Taking d ~ ξ yields δμ0/μ0 ~ 0.05, which is comparable to the smallest relative differences in bulk stiffness that elicit significant changes in cellular motility and differentiation on homogeneous substrates67. This suggests that cellular mechanosensing may operate near the fundamental bounds on precision established in this paper.

An alternative strategy that cells could use to collect additional mechanical information is to perform multiple measurements by actively moving to different locations67,68,69. It is then natural to ask when this approach is preferable to remaining in one place and multiplexing probes with different symmetries. To determine the effectiveness of this strategy, we considered a cell that moves in a straight line and executes an isotropic dipolar probe every body length (Supplementary Note 12). For a 10 μm cell, we find that a cell must move about eight body lengths before δμ0/μ0 becomes smaller than 0.05, the corresponding value for a stationary cell from the previous paragraph. Thus, we expect that a cell would prefer to exert multiple probes in a single location if this could be done in less than the time it takes to move eight body lengths. Indeed, the precision of sensing improves more rapidly with the number of measurements for a stationary sensor than for a moving sensor (Supplementary Note 12). This suggests that cells should choose to stay put and probe their environment to the full extent permitted by their resolution d, after which they have no choice but to begin moving.

The precision of active microrheology

In this section, we apply our modeling framework to active microrheology in disordered polymeric systems and other soft materials. Active microrheology can be used for medical diagnostics, such as monitoring the progression of a cancer52, and also for quality control during material fabrication53,54. Although previous studies have shown that active microrheology can be used to probe the properties of small systems47, it has remained unclear what this technique can learn about heterogeneous materials. What is the fractional uncertainty of active microrheology, and what is the best strategy to estimate the average stiffness of a material using this approach?

To be concrete, we consider active microrheology that harnesses an engineered device to apply a static force and measure displacement at a single location inside an elastic medium (corresponding to the low-frequency limit of a general, frequency-dependent microrheology measurement). In contrast to our considerations regarding cells, such a device may exert a net force on the medium, e.g., by manipulating a bead using magnetic tweezers47. Thus, to explore what such a device can learn beyond what a cell is capable of discerning, we focus on a measurement protocol that consists of a monopole stimulus field and a monopole weight field applied to an elastic network. Such probe fields induce diverging deformations at the points of application. We account for these unphysical divergences by taking the measurement protocol to include a spherical cutoff of radius ξ equal to the mesh size of the network (Supplementary Note 13). For the reconstituted collagen network we considered in the previous section, we used numerical integration to compute the fractional uncertainty of this measurement protocol and found δμ0/μ0 ~ 0.2. This value is small enough to reliably identify the presence of a cancer70, which suggests that even one-particle active microrheology could be an effective tool for medical diagnostics.

Conventionally, active microrheology is done by probing the response of a material in a single direction. In a homogeneous material, this method yields results that do not depend on the direction of the probe. However, an anisotropic probe applied to a heterogeneous material generically yields a response that varies with the direction of the probe. In this case, a more precise measurement could be obtained by probing in different directions and combining the results via sensory multiplexing. Intuitively, the amount of information obtained by these probes may be maximized by spreading out their directions as much as possible. We found that δμ0/μ0 does indeed decrease with the number of samples taken, but with negligible improvements beyond δμ0/μ0 ~ 0.15 (Supplementary Note 13). The precision saturates to a non-zero value because monopoles couple to a fixed spatial extent of the medium regardless of orientation. Thus, the best strategy to probe μ0 by one-particle microrheology is to perform three probes in orthogonal directions, which allows a ~25% improvement in precision over a single monopole probe.

Discussion

An understanding of fundamental bounds on sensing precision and information transmission has a long history of spurring advances in the sciences and engineering, ranging from improvements in telephony driven by Shannon’s initial formulation of information theory to investigations of cell signaling and chemotaxis growing out of the Berg-Purcell limit on concentration sensing40,41,42,71. Studies of the limits of sensor performance are valuable both because they imply design constraints that engineered and evolved systems must satisfy and because determining the optimal performance often uncovers strategies to reach this optimum that can be used to improve performance even if the optimum cannot be attained.

Here, in the spirit of these earlier studies, we have quantified what a physical sensor can learn by probing a heterogeneous material. For media with long-ranged response functions, the smallest possible fractional uncertainty in estimating an average material constant is \(\delta \lambda /{\lambda }_{0} \sim ({\Delta }_{\lambda }^{1/2}/{\lambda }_{0}){(d/a)}^{D/2}{(\xi /a)}^{D/2}\) for a ξ d, \({\lambda }_{0}\gg {\Delta }_{\lambda }^{1/2}\), and D > 1. Remarkably, this relation implies that a finite-sized sensor applied to a standard elastic medium can achieve arbitrarily high precision—in effect, averaging the material constant field λ(r) over an arbitrarily large volume—provided that it can perform multiple measurements down to small enough scales d. This “sensory multiplexing” provides a novel design principle for engineering high precision sensors that would be well-suited for applications on the microscopic scale16,17,18,19,20,21,22,23,72.

In practical terms, our results imply that material properties can be estimated most precisely by making several measurements that each impose different multipole symmetries. Importantly, one does not have to do a very large number of measurements to benefit from this strategy, as simply employing dipole and quadrupole probes can yield substantial improvements. These conclusions have implications both for the design of engineered sensors and for the behavior of living cells. In particular, cells can obtain additional information about the stiffness of their environment by applying multiple probes using forces that vary in space and time. Interestingly, such regulation of probes has been found to influence cellular differentiation73, and cells in culture have recently been observed to vary the spatial symmetries of forces they exert65.

For simplicity, we focused primarily on spherical sensors embedded inside a medium. However, our framework can also be used to study different sensory geometries and motile sensors. Many sensors operate on the boundary of media, including cells grown on flat surfaces27. Moreover, cells in connective tissue can become highly elongated66 and undergo directed migration74, both of which may serve as strategies for overcoming spatial correlations. Future experimental studies will be important to investigate the tradeoffs that cells employ between migrating and varying the angular distribution of their focal adhesions.

Throughout the main text, we assumed a vanishingly small material correlation length ξ, which holds provided that d ξ. Our approach can be readily extended to account for a finite correlation length ξ (Supplementary Note 14). Moreover, although we have focused mostly on a simple scalar version of elasticity, we expect our scaling results to hold for a broad range of media with long-ranged response functions, including the elastic medium in section “The precision of biomechanical sensing”. For short-ranged response functions, the smallest possible fractional uncertainty is \(\delta {\lambda }_{0}/{\lambda }_{0} \sim ({\Delta }_{\lambda }^{1/2}/{\lambda }_{0}){(\xi /a)}^{D/2}\) for a ξ and \({\lambda }_{0}\gg {\Delta }_{\lambda }^{1/2}\), which can be interpreted as the familiar \(1/\sqrt{N}\) scaling of measurement uncertainty for N independent samples, where N ~ (a/ξ)D corresponds to the number of effectively independent subvolumes probed by the sensor (Supplementary Note 14). Finally, we assumed that the elastic properties of the medium within the sensing volume are not significantly mismatched from those of the exterior. Extending our model to account for more complicated constitutive relations and other distributions of the disorder are important directions for future research.

We have focused on athermal materials. For thermal materials, the quantities measured by the sensor fluctuate in time. These fluctuations provide an additional source of temporal noise to the inference process, as well as additional response configurations that can be observed by the sensor. Generalizing our approach to account for these effects would provide a comprehensive physical limit to sensing the properties of materials.

In summary, we have elucidated the perception of material properties in physical space. On small scales, structural heterogeneities place limits on the precision of sensing. We modeled these limits for biopolymer networks and found that they are comparable to the bounds observed for cells in experiment67. Going forward, our theory will guide the design of the next generation of sensors that will be capable of probing materials at the fundamental limits of spatial resolution.

Methods

Numerical minimization of δλ0/λ0

To determine the optimal measurement protocol for a sensor that can apply arbitrary probe fields on its boundary, we used Mathematica’s NMinimize function to search for a global minimum of δλ0/λ0. We performed this minimization over the coefficients \({B}_{k}^{(f)}\) and \({B}_{k}^{(w)}\) using the built-in Nelder–Mead method. The accuracy and precision goals were both chosen to be ϵ = 8, and we took the maximum number of iterations to be \({N}_{\max }=1000\). To explore different local minima, we introduced stochasticity by repeating the minimization for 25 random initial seeds for each choice of the parameter \({k}_{\max }\). For each value of \({k}_{\max }\), the minimum fractional uncertainty \(\delta {\lambda }_{0,\min }/{\lambda }_{0}\) reported was taken to be the minimum of the values found among the 25 trials.