1 Introduction
The capability of aquatic animals to accurately perceive their environment plays a crucial role in their survival. Many fish species employ specialized organs to obtain visual, olfactory and tactile cues from their environment which often complement each other. Predator detection by fish using visual or olfactory cues (Hara Reference Hara1975; Ladich & Bass Reference Ladich and Bass2003; Valentinčič Reference Valentinčič2004) is crucial for providing early warning, since mechanical disturbances may be imperceptible at large distances. On the other hand, sensory organs specialized for detecting mechanical disturbances (Schwartz Reference Schwartz1974) take precedence when fish operate in deep or turbid waters, where visual and other sensory mechanisms may become ineffective. In these situations, the burden of collecting sensory information falls primarily on the ‘lateral line’ organ in fish (Dijkgraaf Reference Dijkgraaf1963; Kroese & Schellart Reference Kroese and Schellart1992; Coombs, Hastings & Finneran Reference Coombs, Hastings and Finneran1996; Coombs & Netten Reference Coombs and Netten2005; Bleckmann & Zelick Reference Bleckmann and Zelick2009). These organs are comprised of hair-like mechanoreceptors called neuromasts (figure 1), which generate neuronal impulses when deflected by either the flow shear (superficial neuromasts – Engelmann et al. (Reference Engelmann, Hanke, Mogdans and Bleckmann2000)) or non-zero pressure gradients (sub-surface ‘canal’ neuromasts – Bleckmann & Zelick (Reference Bleckmann and Zelick2009)). An array of such sensors allows fish to discern both the direction and speed of disturbances generated in their surrounding flow (Chambers et al. Reference Chambers, Akanyeti, Venturelli, Ježov, Brown, Kruusmaa, Fiorini and Megill2014; Asadnia et al. Reference Asadnia, Kottapalli, Miao, Warkiani and Triantafyllou2015).
Figure 1. (a) The lateral line in juvenile zebrafish, with neuromasts visible as bright dots on the body surface (adapted with permission from Sapède et al. (Reference Sapède, Gompel, Dambly-Chaudière and Ghysen2002)). We observe a high density of neuromasts in the head and the tail, with sparser distribution along the midsection. (b) A schematic representation of the distribution of mechanoreceptors along the fish body. (c) The neuromasts bend in response to flow, which generates a neuronal response from sensory cells located at the base (adapted with permission from Kottapalli et al. (Reference Kottapalli, Asadnia, Miao and Triantafyllou2013)).
These flow sensors are distributed in distinctive patterns on the body, with the canal neuromasts distributed evenly along the midline from head to tail (Ristroph, Liao & Zhang Reference Ristroph, Liao and Zhang2015), and superficial neuromasts found in dense clusters near the head and tail, with a sparser distribution along the midsection (figure 1). The fact that they are not distributed uniformly over the body, as well as differences in distribution among species inhabiting different hydrodynamic environments (Atema et al. Reference Atema, Fay, Popper and Tavolga1988; Engelmann et al. Reference Engelmann, Hanke, Mogdans and Bleckmann2000; Bleckmann & Zelick Reference Bleckmann and Zelick2009), suggest that neuromast distribution may be optimized for characterizing hydrodynamic disturbances.
Experimental studies have demonstrated that a well-functioning lateral line is crucial for a range of routine behaviour, such as schooling (Pitcher, Partridge & Wardle Reference Pitcher, Partridge and Wardle1976; Partridge & Pitcher Reference Partridge and Pitcher1980), predator evasion (Blaxter & Fuiman Reference Blaxter, Fuiman, Coombs, Görner and Münz1989), prey detection/capture (Hoekstra & Janssen Reference Hoekstra and Janssen1985), reproduction (Satou et al. Reference Satou, Takeuchi, Nishii, Tanabe, Kitamura, Okumoto and Iwata1994), rheotaxis (Dijkgraaf Reference Dijkgraaf1963; Kanter & Coombs Reference Kanter and Coombs2003), obstacle avoidance (Hassan Reference Hassan, Coombs, Görner and Münz1989) and station-keeping by countering the effects of unsteady gusts (Sutterlin & Waddy Reference Sutterlin and Waddy1975). Disrupting the normal functioning of the lateral line, either via chemical or mechanical means, hinders a fish’s ability to perform these tasks effectively. Liao (Reference Liao2006) demonstrated that disabling the lateral line system influences a fish’s ability to harness energy from unsteady flows. The sensory system also plays a vital role in ‘hydrodynamic imaging’, where fish devoid of visual cues swim past walls and unknown objects repeatedly to form a hydrodynamic ‘map’ of their surroundings (Hassan Reference Hassan, Coombs, Görner and Münz1989; Coombs & Montgomery Reference Coombs and Montgomery1999; Montgomery, Coombs & Baker Reference Montgomery, Coombs and Baker2001; Coombs & Braun Reference Coombs and Braun2003). Certain species such as the blind cave fish, which have evolved degenerated sight, rely heavily on this technique for navigation, and for inferring the shape and size of unfamiliar objects (von Campenhausen, Riess & Weissert Reference von Campenhausen, Riess and Weissert1981; de Perera Reference de Perera2004; Windsor, Tan & Montgomery Reference Windsor, Tan and Montgomery2008).
The lateral line system has inspired the design of artificial sensory arrays, given their potential to transform underwater navigation of robotic vehicles (Yang et al. Reference Yang, Chen, Engel, Pandya, Chen, Tucker, Coombs, Jones and Liu2006, Reference Yang, Nguyen, Chen, Lockwood, Tucker, Hu, Bleckmann, Liu and Jones2010; Ježov et al. Reference Ježov, Akanyeti, Chambers and Kruusmaa2012; Kottapalli et al. Reference Kottapalli, Asadnia, Miao, Barbastathis and Triantafyllou2012; Kruusmaa et al. Reference Kruusmaa, Fiorini, Megill, de Vittorio, Akanyeti, Visentin, Chambers, El Daou, Fiazza and Ježov2014; Asadnia et al. Reference Asadnia, Kottapalli, Miao, Warkiani and Triantafyllou2015; Strokina et al. Reference Strokina, Kämäräinen, Tuhtan, Fuentes-Pérez and Kruusmaa2016; Triantafyllou, Weymouth & Miao Reference Triantafyllou, Weymouth and Miao2016; Kottapalli et al. Reference Kottapalli, Bora, Sengupta, Miao and Triantafyllou2018; Yen, Sierra & Guo Reference Yen, Sierra and Guo2018). Such mechanoreceptors would be a vital addition to the already available suite of visual and acoustic sensors, with the added advantage of low energy-consumption, since they operate via passive mechanical deformation. These vibration-detecting sensors would be crucial for navigation, detection and tracking in low-light conditions, or in scenarios where the use of onboard lights or sonar is undesirable, either for maintaining stealth, or for minimally intrusive observation of animals. Current prototypes of such artificial sensors are based on arrays of pressure transducers (Fernandez et al. Reference Fernandez, Maertens, Yaul, Dahl, Lang and Triantafyllou2011; Venturelli et al. Reference Venturelli, Akanyeti, Visentin, Ježov, Chambers, Toming, Brown, Kruusmaa, Megill and Fiorini2012; Xu & Mohseni Reference Xu and Mohseni2017) and mechanically deforming hair-like structures (Yang et al. Reference Yang, Chen, Engel, Pandya, Chen, Tucker, Coombs, Jones and Liu2006; Tao & Yu Reference Tao and Yu2012; Abdulsadda & Tan Reference Abdulsadda and Tan2013; Dagamseh et al. Reference Dagamseh, Wiegerink, Lammerink and Krijnen2013; DeVries et al. Reference DeVries, Lagor, Lei, Tan and Paley2015; Triantafyllou et al. Reference Triantafyllou, Weymouth and Miao2016).
The importance of the lateral line as an essential sensory organ in fish, and its immense potential for driving the bio-inspired design of artificial sensors, has stimulated numerous experimental and model-based studies. The structure and function of these sensory arrays has been investigated via biological experiments, to characterize their response to pressure differences and object-induced vibrations in water (Gray Reference Gray1984; Denton & Gray Reference Denton and Gray1988; Kroese & Schellart Reference Kroese and Schellart1992; Coombs et al. Reference Coombs, Hastings and Finneran1996; Ćurčić-Blake & van Netten Reference Ćurčić-Blake and van Netten2006). Experiments using artificial fish models have tried to emulate these biological studies, using pressure transducers and hair-like sensors to characterize the frequency and range of oscillating spheres (Montgomery & Coombs Reference Montgomery and Coombs1998) and Kármán vortex streets (Venturelli et al. Reference Venturelli, Akanyeti, Visentin, Ježov, Chambers, Toming, Brown, Kruusmaa, Megill and Fiorini2012). Moreover, there have been a number of mathematical model-based studies, that have combined potential-flow solutions with simplified representations of fish-swimming to study the functioning of the lateral line (Hassan Reference Hassan1992; Franosch et al. Reference Franosch, Hagedorn, Goulet, Engelmann and van Hemmen2009; Bouffanais, Weymouth & Yue Reference Bouffanais, Weymouth and Yue2011; Ren & Mohseni Reference Ren and Mohseni2012; Colvert & Kanso Reference Colvert and Kanso2016). A few of these studies have attempted to infer the optimal arrangement of sensors on rigid objects exposed to various flow conditions. Colvert & Kanso (Reference Colvert and Kanso2016) determined the optimal placement of a single sensor-pair on an elliptical body, moving at different orientations in uniform flow. Ahrari et al. (Reference Ahrari, Lei, Sharif, Deb and Tan2017) used simplified analytical representations to determine optimal sensor-arrangement and -orientation on a rigid hydrofoil, which could best characterize a dipole source with six degrees of freedom in three dimensions.
While model-based studies provide important insight regarding sensing, they suffer from certain drawbacks owing to simplified hydrodynamics, and simplistic representations of fish-swimming (e.g. ellipses and rigid airfoils). Neglecting the effects of viscosity in potential-flow based studies is a notable disadvantage, especially when considering larvae swimming at relatively low Reynolds numbers (
$Re$). Moreover, viscous effects play a substantial role in the operation of the lateral line (Triantafyllou et al. Reference Triantafyllou, Weymouth and Miao2016), given that superficial neuromasts are immersed in the fish’s boundary layer, and canal neuromasts encounter low 
$Re$ flow inside constricted channels. The Reynolds number that animals operate at can also have a considerable impact on the functioning of the lateral line (Webb Reference Webb2014), which cannot be accounted for via inviscid assumptions. The importance of viscous effects has also been demonstrated by Rapo et al. (Reference Rapo, Jiang, Grosenbaugh and Coombs2009), who studied the impact of an oscillating sphere on the boundary layer of a vibrating flat plate, albeit using analytical simplifications to circumvent the high computational cost of three-dimensional numerical simulations. Recent studies using two-dimensional viscous computations have attempted to classify wake patterns behind an oscillating airfoil using artificial neural networks (Alsalman, Colvert & Kanso Reference Alsalman, Colvert and Kanso2018; Colvert, Alsalman & Kanso Reference Colvert, Alsalman and Kanso2018). Using flow sensors placed in the wake of the airfoil, they determine that both the spatial distribution of the sensors as well as the flow variable being measured influence the accuracy for predicting wake characteristics. Here, we investigate the role of hydrodynamics in determining the sensor distribution observed in fish, using two-dimensional Navier–Stokes simulations of self-propelled swimmers to overcome the limitations mentioned above. We determine the optimal spatial distribution of sensors via Bayesian optimal experimental design, and we find that the resulting patterns are closely related to sensory layouts found in natural swimmers.
2 Methods
The present study relies on two-dimensional simulations of a self-propelled swimmer possessing shear stress and pressure gradient sensors on its surface. The swimmer is exposed to disturbances generated by cylinders located at various positions in the environment. The sensor locations are identified by formulating a Bayesian optimal experimental design with the goal of maximizing the information gain of the swimmer in its environment.
2.1 Numerical methods
We conduct two-dimensional simulations of viscous flows past multiple bodies by discretizing the vorticity form of the incompressible Navier–Stokes equations
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D714}=\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D714}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D712}(\boldsymbol{u}_{s}-\boldsymbol{u})),\end{eqnarray}$$ where 
$\boldsymbol{u}$ is the flow-velocity and 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ is the vorticity. The penalty term, 
$\unicode[STIX]{x1D706}\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D712}(\boldsymbol{u}_{s}-\boldsymbol{u}))$ models the interaction of objects with the surrounding fluid (Coquerelle & Cottet Reference Coquerelle and Cottet2008), where 
$0<\unicode[STIX]{x1D712}\leqslant 1$ indicates the solid body. Here 
$\unicode[STIX]{x1D706}$ is the penalization parameter and 
$\boldsymbol{u}_{s}$ represents the combined translational, rotational and deformational velocity of the solid object. The equations are discretized using remeshed vortex methods (Koumoutsakos & Leonard Reference Koumoutsakos and Leonard1995) and wavelet adapted grids (Rossinelli et al. Reference Rossinelli, Hejazialhosseini, van Rees, Gazzola, Bergdorf and Koumoutsakos2015), and the penalty term is integrated via the fully implicit backward Euler method. Additional details for the computational methods may be found in Gazzola et al. (Reference Gazzola, Chatelain, van Rees and Koumoutsakos2011) and Rossinelli et al. (Reference Rossinelli, Hejazialhosseini, van Rees, Gazzola, Bergdorf and Koumoutsakos2015). The simulation domain is a unit square, with an effective resolution of 
$4096^{2}$ grid points. The fish length is 
$L=0.2$ units, with approximately 800 grid points along its midline.
2.2 Swimmer shape and kinematics
We consider two distinct scenarios for the swimmer behaviour to identify the optimal distribution of sensors: one where external disturbances are detected by a static fish-shaped body, and the other involving a self-propelled swimmer. Furthermore, we examine the influence of body geometry on optimal sensor distribution by considering two shapes for the swimmers modelled after zebrafish in their larval and adult stages. The larva shape, shown in figure 2(a), is based on silhouettes extracted from experiments, whereas the adult fish is modelled using a geometric combination of circular arcs, lines and parabolic sections (figure 2b) (Gazzola, Van Rees & Koumoutsakos Reference Gazzola, Van Rees and Koumoutsakos2012). Details regarding shape parameterization for both cases are provided in appendix A: (A 1) and (A 2). The swimmers propel themselves by imposing a sinusoidal wave travelling along the body. Details of the swimming kinematics are also provided in appendix A.
Figure 2. (a) A larva-shaped swimmer detecting disturbances generated by a rotating cylinder (angular velocity 
$=$ 
$10~\text{rotations}~\text{s}^{-1}$). Regions with positive vorticity are coloured in red, and those with negative vorticity are coloured in blue. (b) An adult-shaped swimmer detecting an oscillating cylinder (amplitude 
$=$ 
$0.075L$, frequency 
$=$ 10 Hz). Animations for these two cases are shown in supplementary Movie 1 and Movie 2 and are available online at https://doi.org/10.1017/jfm.2019.940.
2.3 Disturbance-generation and detection
The sensory cues detected by the rigid and swimming bodies described in § 2.2 are generated using oscillating and rotating cylinders of diameter 
$D=0.25L$ (figure 2), and a D-shaped half-cylinder of diameter 
$0.5L$. The amplitude and frequency of the horizontally oscillating cylinders are set to 
$A_{cyl}=0.075L$ and 
$f_{cyl}=10~\text{Hz}$, whereas the angular velocity of the rotating cylinders is set to 
$20\unicode[STIX]{x03C0}~\text{rad}~\text{s}^{-1}$ (
$10~\text{rotations}~\text{s}^{-1}$). The cylinders are placed at various locations within a prescribed region in the computational domain (in the prior-region), as shown in figure 3.
Figure 3. Set-up used for determining the optimal sensor distribution on a fish-like body swimming past a cylinder located within the rectangular area, which is referred to as the prior-region. The sensor placement algorithm attempts to find the best arrangement of sensors that allows the swimmer to identify the correct cylinder position with minimal uncertainty.
We distinguish two types of sensors on the swimmer body. Shear stress sensors estimate the local shear stress by measuring the tangential flow velocity in the reference frame of the swimmer at two grid cells away from the body, corresponding to a physical distance of 
$0.0024L$. These sensors are analogous to superficial neuromasts in fish that protrude into the boundary layer and measure tangential velocity (Kroese & Schellart Reference Kroese and Schellart1992; Bleckmann & Zelick Reference Bleckmann and Zelick2009; Asadnia et al. Reference Asadnia, Kottapalli, Miao, Warkiani and Triantafyllou2015). In addition, we consider pressure gradient sensors that correspond to canal neuromasts observed in natural swimmers. We compute pressure gradient along the swimmers’ surface by fitting a least-squares cubic spline to surface pressure, in order to minimize derivative noise. We note that the sensor measurements are scalar quantities, since both the shear stress and pressure gradient are projected along the surface tangent vector at each measurement location on the body. This is done so that the measured scalar quantities are a close representation of the flow-induced forces that deflect the hair-like sensory structures in real fish.
In the case of self-propelled swimmers, measurements are taken towards the end of a coasting phase to allow self-generated disturbances to subside sufficiently, and are averaged over a small time window from 
$15.750T$ to 
$15.875T$. For motionless larvae, the disturbance sources start moving at 0 s, and time averaging of the recorded data is done between 0.95 s and 1.0 s. This allows transients from the initial cylinder start-up to dissipate sufficiently. Time averaging for the D-cylinder simulations is done from 18 s to 20 s, which allows adequate time for vortex shedding to exhibit a periodically repeating pattern. These measurements are then used to determine the optimal arrangement of sensors on the swimmer body, via the Bayesian optimal experimental design algorithm described in § 2.4.
2.4 Bayesian optimal sensor placement
2.4.1 Bayesian estimation of disturbance location
We consider a disturbance-generating source (for example an oscillating or rotating cylinder) located at coordinates 
$\boldsymbol{r}=(x,y)$ in the region shown in figure 3. The uncertainty in the values of the coordinates of the cylinder is quantified by a probability distribution that is updated based on measurements collected on the surface of the swimmers. The cylinder location can be detected provided that disturbances induced by the cylinder to the surrounding fluid are detected by sensors located on the swimmer surface. The problem of optimal placement implies that we identify the configuration of sensors that can provide the best estimate for the coordinates of the cylinder (
$\boldsymbol{r}$). We assume that the sensor locations are placed symmetrically on both sides of the two-dimensional swimmer and they are described by a vector 
$\boldsymbol{s}\in R^{n}$, that is the midline coordinate of each sensor pair, with values in 
$[0,L]$ (see figure 3). The shear stress or the pressure gradient are measured on the surface points corresponding to the positions 
$\boldsymbol{s}$, and are listed in a vector 
$\boldsymbol{y}\in R^{2n}$.
We denote as 
$\boldsymbol{F}(\boldsymbol{r};\boldsymbol{s})$ the predictions of shear stress/pressure gradient at sensor locations 
$\boldsymbol{s}$, obtained by solving the Navier–Stokes equations with a disturbance-generating source located at 
$\boldsymbol{r}$. Moreover, we assume that we have prior knowledge about the parameter 
$\boldsymbol{r}$, encoded in a prior probability distribution 
$p(\boldsymbol{r})$.
After observing the measurements 
$\boldsymbol{y}$ from sensors 
$\boldsymbol{s}$, we use Bayesian inference to update our prior belief for the plausible values of parameter 
$\boldsymbol{r}$, by identifying the posterior probability distribution 
$p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})$. Following Bayes’ rule, the posterior distribution 
$p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})$ of the model parameters is proportional to the product of the prior distribution 
$p(\boldsymbol{r})$ and the likelihood 
$p(\boldsymbol{y}|\boldsymbol{r},\boldsymbol{s})$. The likelihood function represents the probability that a particular measurement 
$\boldsymbol{y}$ for a given sensor arrangement 
$\boldsymbol{s}$ originates from the disturbance source located at 
$\boldsymbol{r}$. We assume a prediction error, 
$\unicode[STIX]{x1D73A}(\boldsymbol{s})$, as the difference between the measurements 
$\boldsymbol{y}$ and the predictions 
$\boldsymbol{F}(\boldsymbol{r};\boldsymbol{s})$ such that
$$\begin{eqnarray}\boldsymbol{y}=\boldsymbol{F}(\boldsymbol{r};\boldsymbol{s})+\unicode[STIX]{x1D73A}(\boldsymbol{s}).\end{eqnarray}$$ The prediction error term (
$\unicode[STIX]{x1D73A}(\boldsymbol{s})$) represents errors that can be attributed to measurement- and model-errors, as well as numerical errors due to spatio-temporal discretization of the Navier–Stokes equations. Following the maximum entropy criteria the prediction error 
$\unicode[STIX]{x1D73A}(\boldsymbol{s})$ follows a multivariate Gaussian distribution 
${\mathcal{N}}(0,\unicode[STIX]{x1D6F4}(\boldsymbol{s}))$ with zero mean and covariance matrix 
$\unicode[STIX]{x1D6F4}(\boldsymbol{s})\in R^{2n\times 2n}$. The likelihood function 
$p(\boldsymbol{y}|\boldsymbol{r},\boldsymbol{s})$ is then expressed as
$$\begin{eqnarray}\displaystyle p(\boldsymbol{y}|\boldsymbol{r},\boldsymbol{s})=\frac{1}{\sqrt{(2\unicode[STIX]{x03C0})^{n}\det (\unicode[STIX]{x1D72E}(\boldsymbol{s}))}}\exp \left(-\frac{1}{2}(\boldsymbol{y}-\boldsymbol{F}(\boldsymbol{r};\boldsymbol{s}))^{\text{T}}\unicode[STIX]{x1D72E}^{-1}(\boldsymbol{s})(\boldsymbol{y}-\boldsymbol{F}(\boldsymbol{r};\boldsymbol{s}))\right). & & \displaystyle\end{eqnarray}$$2.4.2 Optimal sensor placement based on information gain
The goal of the optimal sensor placement problem is to find the locations 
$\boldsymbol{s}$ of the sensors such that the data measured in these locations are most informative for estimating the position 
$\boldsymbol{r}$ of the disturbance. A measure of information gain is provided by the Kullback–Leibler divergence between the prior and the posterior distribution. We postulate that the optimal sensor configuration maximizes a utility function that represents the information gain, or equivalently, the Kullback–Leibler divergence defined as
$$\begin{eqnarray}u(\boldsymbol{s},\boldsymbol{y}):=\int _{{\mathcal{R}}}p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})\ln \frac{p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})}{p(\boldsymbol{r})}\,\text{d}\boldsymbol{r}.\end{eqnarray}$$ We note that in the experimental design phase, the measurements 
$\boldsymbol{y}$ are not available. Thus, the prediction error model (2.2) is used to generate measurements for given model parameter values 
$\boldsymbol{r}$ and sensor configuration 
$\boldsymbol{s}$. We identify the best sensor arrangement by maximizing a utility function, defined as the expected value of the Kullback–Leibler divergence over all possible values of the measurements simulated by (2.2) (Ryan Reference Ryan2003):
$$\begin{eqnarray}\displaystyle U(\boldsymbol{s}):=\mathbb{E}_{\boldsymbol{y}|\boldsymbol{s}}[u(\boldsymbol{s},\boldsymbol{y})] & = & \displaystyle \int _{{\mathcal{Y}}}u(\boldsymbol{s},\boldsymbol{y})\;p(\boldsymbol{y}|\boldsymbol{s})\,\text{d}\boldsymbol{y}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\mathcal{Y}}}\int _{{\mathcal{R}}}p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})\ln \frac{p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})}{p(\boldsymbol{r})}\;p(\boldsymbol{y}|\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}.\end{eqnarray}$$ The expected utility function involves a double integral over the parameter space 
$\boldsymbol{r}$ and over the measured data 
$\boldsymbol{y}$. An efficient estimator of this double integral using sampling techniques is provided by Huan & Marzouk (Reference Huan and Marzouk2013). A similar estimator is used in the present work,
$$\begin{eqnarray}\hat{U} (\boldsymbol{s})=\frac{1}{N_{\boldsymbol{y}}}\mathop{\sum }_{j=1}^{N_{\boldsymbol{y}}}\mathop{\sum }_{i=1}^{N_{\boldsymbol{r}}}w_{i}p(\boldsymbol{r}^{(i)})\left[\ln p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r}^{(i)},\boldsymbol{s})-\ln \left(\mathop{\sum }_{k=1}^{N_{\boldsymbol{r}}}w_{k}p(\boldsymbol{r}^{(k)})p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r}^{(k)},\boldsymbol{s})\right)\right].\end{eqnarray}$$ A detailed derivation and discussion of the estimator is provided in appendix B. Our estimator employs a quadrature technique to evaluate the integral over the two-dimensional parameter space 
$\boldsymbol{r}$. In (2.6), 
$\boldsymbol{r}^{(i)}$ and 
$w_{i}$ denote 
$N_{\boldsymbol{r}}$ quadrature points and corresponding weights related to discretization of the two-dimensional prior-region (figure 3). A total of 
$N_{\boldsymbol{r}}$ distinct Navier–Stokes simulations are conducted, with a cylinder positioned at various discrete points 
$\boldsymbol{r}^{(i)}$, and the quadrature is evaluated using the trapezoidal rule. Based on the prediction error defined in (2.2), the measured data 
$\boldsymbol{y}^{(i,j)}$ in (2.6) are given by
$$\begin{eqnarray}\boldsymbol{y}^{(i,j)}=\boldsymbol{F}(\boldsymbol{r}^{(i)};\boldsymbol{s})+\unicode[STIX]{x1D73A}^{(j)},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D73A}^{(j)}$, with 
$j=1,\ldots ,N_{\boldsymbol{y}}$, are vectors sampled from the distribution 
${\mathcal{N}}(0,\unicode[STIX]{x1D6F4}(\boldsymbol{s}))$. 
$N_{\boldsymbol{y}}$ is set to 
$100$ in the current work, which results in a smoother estimate of 
$\hat{U} (\boldsymbol{s})$ in (2.6).
We note that the computational effort for evaluating 
$\hat{U} (\boldsymbol{s})$ in (2.6) depends primarily on the number of Navier–Stokes simulations, 
$N_{\boldsymbol{r}}$, which are required to evaluate 
$\boldsymbol{F}(\boldsymbol{r}^{(i)},\boldsymbol{s})$ for different disturbance locations 
$\boldsymbol{r}^{(i)}$, and subsequently to determine 
$\boldsymbol{y}^{(i,j)}$ using (2.7). The computational burden does not depend on the number of measured samples 
$N_{\boldsymbol{y}}$, since there are no additional time-consuming simulations involved in generating 
$\unicode[STIX]{x1D73A}^{(j)}$. Thus the computational effort scales linearly with the number 
$N_{r}$ of model parameter points 
$\boldsymbol{r}^{(i)}$.
We assume that the prior distribution 
$p(\boldsymbol{r})$ for the location of the disturbance source is uniform over the prior-region shown in figure 3, i.e. the probability of finding the source is constant for all locations. Moreover, the only available information we have is a description of the prior-region where the disturbance may be found. Using Bayes’ theorem, and the fact that the prior distribution is uniform, we can assert that the posterior distribution of a disturbance location 
$\boldsymbol{r}$, 
$p(\boldsymbol{r}|\boldsymbol{y},\boldsymbol{s})$, is proportional to the likelihood function 
$p(\boldsymbol{y}|\boldsymbol{r},\boldsymbol{s})$.
The covariance matrix 
$\unicode[STIX]{x1D6F4}(\boldsymbol{s})$ depends primarily on the sensor positions 
$\boldsymbol{s}$, and is diagonal if the errors at the given sensor positions are independent of each other. In the current work, the prediction errors are assumed to be correlated for measurements collected on the same side of the swimmer (i.e. left- or right-lateral surfaces), and decorrelated if the measurements originate from opposite sides. An exponentially decaying correlation is assumed for the covariance matrix
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{ij}(\boldsymbol{s})=\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D70E}^{2}\exp \left(-\frac{\Vert \boldsymbol{x}(s_{i})-\boldsymbol{x}(s_{j})\Vert }{\ell }\right),\quad & \text{if }1\leqslant i,j\leqslant n,\\ \unicode[STIX]{x1D6F4}_{i-n,j-n}(\boldsymbol{s}),\quad & \text{if }n<i,j\leqslant 2n,\\ 0,\quad & \text{otherwise,}\end{array}\right.\end{eqnarray}$$ where 
$\boldsymbol{x}(s_{i})$ corresponds to the coordinates of the 
$i$th sensor on the right-lateral surface of the swimmer, 
$\ell >0$ is the prescribed correlation length and 
$\unicode[STIX]{x1D70E}$ is the correlation strength. For all the simulations described in this work, the correlation length is set to be 
$\ell =0.01L$. The correlation strength 
$\unicode[STIX]{x1D70E}$ is a fixed percentage (
$30\,\%$) of the mean sensor-measurement, which is computed over all available instances of 
$\boldsymbol{r}$ and at all points discretizing the swimmer skin. This form of the correlation error reduces the information gain when sensors are placed too close together (Papadimitriou & Lombaert Reference Papadimitriou and Lombaert2012; Simoen, Papadimitriou & Lombaert Reference Simoen, Papadimitriou and Lombaert2013), and prevents excessive clustering of sensors within confined neighbourhoods.
Finally, we provide an intuitive interpretation of how (2.6) relates to information gain. Let us assume that a particular set of sensors is able to characterize the disturbance sources quite effectively. Moreover, we assume that the measurement 
$\boldsymbol{y}^{(i,j)}$ has been generated by a disturbance located at 
$\boldsymbol{r}^{(i)}$. This implies that the posterior 
$p(\boldsymbol{r}|\boldsymbol{y}^{(i,j)},\boldsymbol{s})$, which indicates the probability that a particular disturbance source 
$\boldsymbol{r}$ has generated the measurement 
$\boldsymbol{y}^{(i,j)}$, is peaked and centred around the true source location 
$\boldsymbol{r}^{(i)}$. Since the prior distribution is uniform, the likelihood 
$p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r},\boldsymbol{s})$ is proportional to the posterior, and is also peaked and centred around 
$\boldsymbol{r}^{(i)}$. Thus, the first term in (2.6) is large, whereas most of the terms in the second sum are close to zero, since the probability of measurement 
$\boldsymbol{y}^{(i,j)}$ originating from source 
$\boldsymbol{r}^{(k)}$ is small due to the peaked nature of the posterior (except for 
$k=i$). In this case, the expected utility value computed using (2.6) is large. On the other hand, a poor sensor arrangement which cannot characterize source positions well, yields flatter likelihood and posterior distributions due to high uncertainty. Thus, different source positions yield similar measurements at the selected sensors, which makes the second sum in (2.6) larger (non-zero 
$p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r}^{(k)},\boldsymbol{s})$ even for 
$k\neq i$), thereby reducing the utility value.
2.4.3 Optimization of the expected utility function
The optimal sensor arrangement is obtained by maximizing the expected utility estimator 
$\hat{U} (\boldsymbol{s})$ described in (2.6). However, optimal sensor placement problems are characterized by a relatively large number of multiple local optima. Heuristic approaches, such as the sequential sensor placement algorithm described by Papadimitriou (Reference Papadimitriou2004), have been demonstrated to be effective alternatives. In this approach, the optimization is carried out iteratively, one sensor at a time. First, 
$\hat{U} (\boldsymbol{s})$ is computed for a single sensor-pair 
$\boldsymbol{s}=s_{1}$, and the optimal solution 
$s_{1}^{\star }$ is obtained by identifying the maximum in 
$\hat{U} (\boldsymbol{s})$. Then, 
$\hat{U} (\boldsymbol{s})$ is recomputed with 
$\boldsymbol{s}=(s_{1}^{\star },s_{2})$, and it is optimized with respect to the second sensor-pair, resulting in an optimal solution 
$s_{2}^{\star }$. We can generalize this procedure for all subsequent sensors, by defining 
$\hat{U} _{i}(s)=\hat{U} (s_{1}^{\star },\ldots ,s_{i-1}^{\star },s)$. The optimal solution for the 
$i$th sensor is given as
$$\begin{eqnarray}s_{i}^{\star }=\text{argmax}_{s}\hat{U} _{i}(s)\quad \text{and}\quad \hat{U} _{i}^{\star }=\max _{s}\hat{U} _{i}(s).\end{eqnarray}$$ We note that the scalar variable 
$s$ denotes the position of a single sensor-pair, whereas the vector 
$\boldsymbol{s}$ holds the position of all sensor-pairs along the swimmer’s midline.
The sequential placement procedure is carried out for a number of sensors, 
$N_{s}$, and it terminates when the last sensor in the optimal configuration is identified 
$\boldsymbol{s}^{\star }=(s_{1}^{\star },\ldots ,s_{N_{s}}^{\star })$. Papadimitriou (Reference Papadimitriou2004) has demonstrated that the heuristic sequential sensor placement algorithm provides a sufficiently accurate approximation of the global optimum. Moreover, using the sequential optimization approach, 
$N_{s}$ one-dimensional problems have to be solved, instead of one 
$N_{s}$-dimensional problem. We solve each one-dimensional problem of identifying the maximum of 
$\hat{U} _{i}$ via a grid search, where the swimmer midline is discretized using the points 
$\{k\unicode[STIX]{x0394}s,k=0,\ldots ,N_{g}\}$, with 
$\unicode[STIX]{x0394}s=L/N_{g}$ and 
$N_{g}=1000$. Thus, for each iteration of sequential optimization, the utility estimator in (2.6) has to be evaluated 
$N_{g}+1$ times.
We remark that the Bayesian optimal design procedure is computationally demanding, as it entails model simulations for several different sensor configurations 
$\boldsymbol{s}$. To minimize the relevant computational cost, we run 
$N_{\boldsymbol{r}}$ distinct Navier–Stokes simulations for all disturbance locations 
$\boldsymbol{r}^{(i)}$ (
$i=1\ldots ,N_{r}$), and store the shear stress and pressure gradient at all available discretization points along the swimmer skin offline. This allows us to reuse simulation data for a particular disturbance source, without having to re-run Navier–Stokes simulations for different sensor configurations. We note that the skin discretization may not correspond to the 
$N_{g}$ points used for computing 
$\hat{U} _{i}^{\star }$. Thus, the output quantities of interest are averaged at appropriate locations along the swimmer surface, over a small neighbourhood of size 
$0.01L$.
Figure 4. Snapshots of the vorticity field around a static larva profile in the presence of a horizontally oscillating cylinder. The snapshots are taken at regular intervals over a single oscillation period, with positive vorticity shown in red and negative vorticity shown in blue. A corresponding animation is shown in supplementary movie 3.
Figure 5. Utility plots for a stationary, larva-shaped body with (a) oscillating and (b) rotating cylinders. The curves indicate the utility for placing the first shear stress sensor at a given location 
$s$. The utility curves were not computed in the region 
$0.95<s/L\leqslant 1$, to avoid potential numerical issues resulting from sharp corners at the tail. (c,e) Standard deviation of horizontal and vertical velocity caused by oscillating cylinders, with larger deviation shown in yellow and lower values shown in black. The standard deviation was computed across nine distinct simulations (six time-snapshots recorded in each simulation), with a single oscillating cylinder placed at nine locations uniformly in the prior-region. (d,f) Standard deviation of velocity components for the rotating cylinders. (g) Optimal sensor distribution determined using sequential placement. Sensors for detecting oscillating cylinders are shown as black squares, whereas those for detecting rotating cylinders are shown as blue circles. The numbering indicates the sequence determined by the optimal placement algorithm.
3 Results
We first examine the optimal arrangement of shear stress and pressure gradient sensors on motionless larva in the presence of oscillating, rotating and D-shaped cylinders. We then consider self-propelled swimmers, which are exposed to cylinder-generated disturbances.
3.1 Stationary swimmer in the vicinity of oscillating/rotating cylinders
We first consider the set-up of a stationary larva-shaped swimmer and a cylinder that either oscillates parallel to the ‘anteroposterior’ axis of the body, or rotates with a constant angular velocity. The oscillating-cylinder set-up is shown in figure 4, and depicts the vorticity generated by the cylinder along the larva’s body. These two set-ups allow us to analyse mechanical cues (i.e. vibrations in the flow field) without interference from a self-generated boundary layer, which has a tendency to obscure external signals in the case of towed and self-propelled bodies. The simulation domain extends from 
$[0,1]$ in 
$x$ and 
$y$, and the rectangular prior-region in both the set-ups corresponds to 
$\boldsymbol{r}_{min}=(0.357,0.375)$ and 
$\boldsymbol{r}_{max}=(0.7,0.47)$. A total of 
$11\times 37=407$ potential 
$\boldsymbol{r}^{(i)}$ locations are distributed uniformly throughout the region, and the static object’s centre of mass is located at 
$(0.5,0.3)$. The kinematic viscosity is set to 
$\unicode[STIX]{x1D708}=1\times 10^{-4}$ in these simulations.
3.1.1 The utility function, and sensor placement
The optimal distribution of sensors along the larva’s body can be determined using the estimator 
$\hat{U} (\boldsymbol{s})$ defined in (2.6). Higher utility values indicate that measurements taken at the corresponding locations are more informative. More specifically, the utility at location 
$s$ is high if a sensor placed there can more effectively differentiate between signals originating from distinct cylinder locations. The utility curves computed from signals generated by the oscillating and rotating cylinders are shown in figures 5(a) and 5(b). The maxima in these curves suggest that the best location for detecting both oscillating and rotating cylinders, using shear stress sensors, is at 
$s/L=0.033$. We remark that this location involves a notable change in body-surface curvature, as can be discerned from the swimmer silhouettes shown in the figures.
We postulate that the best sensor positions are those that are exposed to large variations in the quantity of interest, namely the shear stress or pressure gradient, since this would allow the sensors to best distinguish between different disturbance sources more readily. We confirm that this is indeed the case, by visualizing the standard deviation of velocity components in regions surrounding the larva, in figures 5(c)–5(f). The standard deviation measures the variation among simulations when cylinders are placed at different positions in the prior-region shown in figure 3. The colour scales are identical for panels 5(c) and 5(e) (the oscillating cylinder scenario), but different from the colour scales in panels 5(d) and 5(f) (the rotating cylinder scenario). The colour scales in 5(d) and 5(f) have been saturated by approximately 30 times, so that weaker flow disturbances created by the rotating cylinders are adequately visible.
We observe from figures 5(c) and 5(e) that changing the position of an oscillating cylinder gives rise to significant differences in the tangential velocity (shear stress) close to the head and the tail. This implies that signals measured by sensors in these regions differ markedly from one simulation to the other, which arguably would make it easier to estimate the position of a particular cylinder. A large variation in horizontal velocity 
$u$ occurs close to a change in body curvature at 
$s/L\approx 0.033$, which also corresponds to the global maximum in 
$\hat{U} _{1}(s)$ (figure 5a). The utility curve exhibits consistently high values for 
$s/L\leqslant 0.15$, which results from large variations in 
$u$ and 
$v$ in regions surrounding the head. We note that large variations in the lateral velocity 
$v$ occur primarily at the head- and tail-tip (figure 5e), with almost no variation along the midsection (
$0.2<s/L\leqslant 1$). This can be attributed to 
$v$ being almost zero in these regions (across all simulations), owing to negligible recirculation along these relatively straight body sections. The large variation in 
$v$ at the head/tail tip may be explained by the flow turning at the corners, as is evident from the time-series snapshots shown in figure 4. We note that while 
$u$ appears to exhibit large deviation around the midsection (
$0.4\leqslant s/L\leqslant 0.6$ in figure 5c), the utility curve in figure 5(a) does not show a corresponding spike. This may be related to the fact that the standard deviation plots were compiled using a small subset of nine cylinder locations out of the 407 used for the utility plot. Furthermore, a close inspection of figure 5(c) indicates that these large deviations in 
$u$ near the midsection occur beyond the detection range of the sensors, i.e. too far away to be picked up by microscopic neuromasts that are 
$0.0024L$ in length.
As in the case of oscillating cylinders, the standard deviation plots for rotating cylinders in figures 5(d) and 5(f) can be correlated to the utility curve in figure 5(b); high utility values (
$s/L\leqslant 0.15$, figure 5b) correspond to large deviations in both 
$u$ and 
$v$ near the head (figure 5d,f). Based on the utility curve, the highest sensitivity for measuring flow perturbations corresponds to the head and posterior sections of the body. This suggests that the head and tail are the most informative regions for detecting shear stress fluctuations for a static larva, regardless of the type of disturbance being considered. This observation is consistent with the distribution of neuromasts shown in figure 1, where the surface neuromasts are visible in high concentrations in the head and posterior regions of fish, but show sparse presence along the midsection.
Figure 6. (a) Utility curves for placing the first three sensors on a static larva that detects oscillating cylinders (figure 4). The solid green curve corresponds to 
$\hat{U} _{1}(s)$, the dashed purple curve to 
$\hat{U} _{2}(s)$ and the red dash-dot curve to 
$\hat{U} _{3}(s)$. (b) The optimal utility 
$\hat{U} _{n}^{\star }$ for the 
$n$th sensor can be determined as 
$\max _{s}\hat{U} _{n}(s)$ using the curves shown in panel (a) (see also (2.9)).
3.1.2 Sequential sensor placement
In the previous section we discussed the case of a single sensor on the swimmer body. We now examine the optimal arrangement of multiple sensors, where the best location for the 
$n$th sensor is determined provided that 
$n-1$ sensors have already been placed. Assume that the first sensor has been placed at 
$s_{1}^{\star }$ using the global maximum in utility curve 
$\hat{U} _{1}(s)$. The next best sensor-location is determined by recomputing the utility function 
$\hat{U} _{2}(s)$ as described in § 2.4.3. Following this procedure, the optimal location of all sensors is determined sequentially.
Figure 5(g) shows the optimal distribution of 20 sensors for the static larva determined in this manner. We first examine the optimal arrangement for detecting oscillating cylinders, with the corresponding sensors depicted as black squares. We observe that out of the first 10 sensors, numbers 1, 3, 5, 9 are placed at the head, whereas numbers 2, 4, 6, 7, 8, 10 are found towards the posterior. This suggests a large information gain via sensors located in the head and the tail. For detecting rotating cylinders (sensors shown as blue circles in figure 5g), sensors 1, 3, 4, 6, 8, 10 are found in the head, and sensors 2, 5, 7, 9 are placed in the posterior section.
We also examine the utility curves for placing the first three oscillation-detecting shear stress sensors in figure 6(a). We observe that 
$\hat{U} _{2}(s_{1})\approx \hat{U} _{1}(s_{1})$, which indicates that placing a second sensor at the same location as the first (
$s_{1}/L=0.033$) would not lead to an appreciable increase in the utility value (i.e. no gain in useful information). The maximum in 
$\hat{U} _{2}(s)$ occurs at 
$s/L=0.95$, which yields the optimal location 
$s_{2}^{\star }$ for the second sensor. Another notable aspect of curve 
$\hat{U} _{2}(s)$ is a pronounced ‘v-shaped’ depression in the vicinity of 
$s_{1}^{\star }$, which results from using a non-zero correlation length in (2.8). The low utility values in this region impede the placement of sensors too close to each other. Using a zero correlation length would have resulted in an abrupt drop in 
$\hat{U} _{2}(s)$ at 
$s_{1}^{\star }$ (instead of the smooth depression), and could lead to excessive clustering of sensors within confined neighbourhoods. Figure 6(b) shows the cumulative utility value for an increasing number of sensors placed on the swimmer body. We observe that after a rapid initial rise for the first three to five sensors, the utility of placing subsequent sensors increases very slowly. This indicates that using a limited number of optimal sensor locations should be sufficient to characterize disturbance sources with reasonably good accuracy.
3.2 Motionless larva in the wake of a D-cylinder
We now consider simulations where a rigid larva-shaped profile is placed in the unsteady vortex-wake generated by a D-shaped half-cylinder (figure 7). This configuration is inspired by the pioneering work of Liao et al. (Reference Liao, Beal, Lauder and Triantafyllou2003) who examined the fluid dynamics of trout placing themselves behind rocks. A uniform horizontal flow of 
$1L~\text{s}^{-1}$ is imposed throughout the computational domain, and the rigid bodies are held stationary. The D-cylinder is located at 
$(0.2,0.5)$, and the rectangular prior-region for placing the larvae extends from 
$\boldsymbol{r}_{min}=(0.3,0.43)$ to 
$\boldsymbol{r}_{max}=(0.79,0.57)$. A total of 
$11\times 36=396$ potential 
$\boldsymbol{r}^{(i)}$ locations are distributed uniformly throughout the prior-region. The Reynolds number is 
$Re=200$ based on the cylinder diameter, and 
$Re=400$ based on the swimmer length.
Figure 7. Snapshots of the vorticity field around a static larva in the wake of a D-shaped cylinder with diameter 
$0.5L$. The D-cylinder is oriented at a 
$10^{\circ }$ angle with respect to a uniform horizontal flow to promote vortex shedding. The snapshots are shown at regular time intervals, with positive vorticity shown in red and negative vorticity shown in blue. A corresponding animation is shown in Movie 4.
Figure 8 shows the utility curve for placing the first shear stress sensor on the static larva, as well as the sensor distribution resulting from sequential placement. The utility values for 
$0.2\leqslant s/L\leqslant 0.6$ are close to zero, which implies that placing the first sensor along the midsection would provide minimal information gain. Using the sequential-placement procedure described in § 2.4.3, we determine that all of the first 10 sensors are placed at the head, with no sensors present in the tail. Our results indicate that sensors at the head are far more significant than sensors in the mid- and posterior-sections of the body for detecting the unsteady wake behind a half-cylinder.
Figure 8. (a) Utility curve 
$\hat{U} _{1}(s)$ for larvae in a D-cylinder’s wake. (b) Sequential placement of 20 sensors, with the order of placement shown.
3.3 Self-propelled swimmers: shear stress sensors
Fish generate vorticity on their bodies by their undulatory motion. Their flow-sensing neuromasts are completely immersed in this self-generated flow field, which likely has a significant impact on their ability to detect external disturbances. To include the influence of these self-generated flows on optimal sensor placement, we now consider simulations of self-propelled swimmers that are exposed to oscillating and rotating cylinders (figure 2). These swimmers utilize an intermittent swimming gait referred to as ‘burst and coast’ swimming, which allows for improved sensory perception (Kramer & McLaughlin Reference Kramer and McLaughlin2001), as self-generated disturbances subside during the coasting phase. The swimmers perform four full burst-coast swimming cycles starting from rest, before the cylinder starts oscillating or rotating, as depicted in movie 1 and movie 2. In the initial transient phase, the swimmer gains a speed of approximately 
$0.7L~\text{s}^{-1}$, which corresponds to a Reynolds number of 
$\mathit{Re}=uL/\unicode[STIX]{x1D708}\approx 280$ (with 
$L=0.2$ and 
$\unicode[STIX]{x1D708}=1\times 10^{-4}$). At the start of the fifth coasting phase, the cylinder starts moving, which simulates the startle/attack response of a prey/predator present in the swimmer’s vicinity. The rectangular prior-region for initializing the cylinders extends from 
$\boldsymbol{r}_{min}=(0.25,0.375)$ to 
$\boldsymbol{r}_{max}=(0.7,0.5)$, with a total of 
$11\times 37=407$ potential 
$\boldsymbol{r}^{(i)}$ locations distributed uniformly throughout the region. The swimmer’s centre of mass is located at 
$(0.5,0.3)$.
To determine the extent to which body shape influences optimal placement, we perform simulations using a larva-shaped profile, and a simplified model of an adult. Figure 9 compares the utility curves and sensor distributions for these two distinct swimmers. Based on the utility curves in figures 9(a) and 9(c), we deduce that the head is the most suitable region for placing the first sensor, as was the case for the motionless profiles examined in the previous sections. We also observe that the utility curves are correlated to the surface curvature of their respective body profiles; in the case of the larva, there is marked variation in 
$\hat{U} _{1}(s)$ for 
$s/L\leqslant 0.2$, which corresponds to large curvature changes in the body surface. The utility curve also shows a gradual variation for 
$s/L\geqslant 0.6$, which corresponds to a gentler change in curvature of the surface. Similarly, the utility curves and body curvature for the adult vary rapidly for 
$s/L\leqslant 0.05$ and more gradually for 
$s/L\geqslant 0.6$. Furthermore, we note that the blue utility curves in figures 9(a) and 9(c) are close to 
$0$ for 
$s/L\geqslant 0.2$. This suggests that the head is the most useful region for placing rotation-detecting shear stress sensors, irrespective of differences in body shape.
Figure 9. (a) Utility curves for the first shear stress sensor, 
$\hat{U} _{1}(s)$, on a larva-shaped swimmer (black squares, oscillating cylinders; blue circles, rotating cylinders). (b) Sequential placement of 20 sensors along the body, with the order of placement shown. (c) Utility curves for an adult-shaped swimmer. (d) Sensor placement for the adult, with results from horizontal and rotating disturbances shown separately for clarity.
The optimal sensor arrangements observed in figures 9(b) and 9(d) are listed in table 1. There is strong indication that the head is the most important region for detecting shear stress caused by external disturbances, followed by the posterior section; the midsection appears to be insensitive to shear-stress variations altogether, as evidenced by the lack of sensors in this region. Moreover, the posterior section appears to be insensitive to rotating disturbances regardless of the body shape. These observations agree well with surface neuromast distributions observed in live fish (figure 1a), where large numbers are found in the head and the tail, with sparser clustering in the midsection.
Table 1. Optimal distribution of the first 10 shear stress sensors for the self-propelled swimmers. The body has been divided into three distinct segments: the head (
$0\leqslant s/L<0.2$), the midsection (
$0.2\leqslant s/L<0.6$) and the posterior (
$0.6\leqslant s/L\leqslant 1$).
3.4 Optimal sensor placement using combined datasets
Fish are subject to a multitude of external stimuli over the course of their lifetime. Hence, it is conceivable that neuromasts may be attuned to diverse sources of disturbance. We emulate this situation for optimal sensor placement by considering data collected from the five different simulation set-ups simultaneously, namely, motionless larvae with oscillating, rotating and D-shaped cylinders, and self-propelled larvae with oscillating and rotating cylinders. The sequential-placement procedure described in the previous sections is followed, with a slight modification to the definition of the utility function. The combined utility for the five different configurations may be expressed as a sum of the individual utility functions, due to the conditional independence of the measurements on the sensor locations (see appendix C). The resulting utility curve for the shear stress sensors is shown in figure 10, along with the optimal sensor distribution. We observe predominant placement of sensors in the head and tail, corresponding to large utility values in these regions. Moreover, we find that virtually no sensors are located in the midsection. The dense clustering of sensors in the head and tail, with sparse distribution in the midsection, yet again resembles surface neuromast patterns found in live fish (figure 1), and indicates that fish extremities may be ideal for detecting variations in shear stress.
Figure 10. (a) Utility curve for the first shear stress sensor, 
$\hat{U} _{1}(s)$, on a larva-shaped swimmer, using a combination of all five flow configurations described in the paper. (b) Sequential placement of 20 shear stress sensors along the body.
Figure 11. (a) Utility curve for the first pressure gradient sensor, 
$\hat{U} _{1}(s)$, on a larva-shaped swimmer, using a combination of all five flow configurations described in the paper. (b) Sequential placement of 20 pressure gradient sensors along the body.
3.5 Optimal pressure gradient sensors
We now consider the optimal placement of pressure gradient sensors on the larva’s body. These sensors are analogous to canal neuromasts found in live fish, which display markedly similar distribution patterns across a variety of fish species (Ristroph et al. Reference Ristroph, Liao and Zhang2015). The canal is usually present in a continuous line running from head to tail, and shows a high concentration of neuromasts in canal branches found in the head (Coombs, Janssen & Webb Reference Coombs, Janssen and Webb1988; Ristroph et al. Reference Ristroph, Liao and Zhang2015). We use a combination of the five distinct flow configurations, described earlier, to determine the optimal arrangement of pressure gradient sensors by following the procedure described in § 3.4. The resulting utility curve and sensor distribution are shown in figure 11. The most notable difference between the arrangement of pressure gradient sensors (figure 11b), and that of shear stress sensors (figure 10b), is observed in the midsection of the body. We find a consistent distribution of pressure gradient sensors in the midsection, which is not the case for shear stress sensors. Out of the 20 pressure gradient sensors placed, 10 are found clustered densely in the head (
$s/L\leqslant 0.1$, which corresponds to high utility values in figure 11a), and the other 10 are spaced regularly throughout the body. This arrangement is similar to the neuromast distribution found in subsurface canals, which yet again suggests that this sensory structure may have evolved for detecting changes in pressure gradients with high accuracy. In fact, the utility curve shown in figure 11(a) agrees qualitatively with the canal density reported by Ristroph et al. (Reference Ristroph, Liao and Zhang2015), especially for 
$s/L<0.2$. However, a direct comparison must be made with care, given that our simulations are two-dimensional, whereas the distributions reported by Ristroph et al. (Reference Ristroph, Liao and Zhang2015) display significant three-dimensional branching in the head.
3.6 Inference of disturbance-generating source
Having determined the optimal distribution of sensors on the swimmer body, we now assess how effectively these arrangements can characterize the disturbance sources. For a given set of sensors 
$\boldsymbol{s}$, this involves estimating the probability that a particular sensor measurement may originate from different cylinder positions within the prior-region. For this, we consider the measurements 
$\boldsymbol{y}^{(GT)}$ at the sensor locations, generated from a single cylinder located at 
$\boldsymbol{r}^{(GT)}$ (the superscript 
$GT$ denotes ‘ground-truth’). For a given sensor configuration 
$\boldsymbol{s}$, the measurements 
$\boldsymbol{y}^{(GT)}$ are computed using the prediction error model 
$\boldsymbol{y}^{(GT)}=\boldsymbol{F}(\boldsymbol{r}^{(GT)};\boldsymbol{s})+\unicode[STIX]{x1D73A}(\boldsymbol{s})$, where 
$\boldsymbol{F}(\boldsymbol{r}^{(GT)};\boldsymbol{s})$ is obtained by simulating the Navier–Stokes equations with an oscillating cylinder located at 
$\boldsymbol{r}^{(GT)}$, and 
$\unicode[STIX]{x1D73A}(\boldsymbol{s})$ is a vector sampled from the Gaussian distribution 
${\mathcal{N}}(0,\unicode[STIX]{x1D6F4}(\boldsymbol{s}))$. Assuming that the disturbance position 
$\boldsymbol{r}^{(GT)}$ is unknown, the swimmer attempts to identify it by assigning probability values to all possible cylinder locations 
$\boldsymbol{r}$ within the prior-region (i.e. by determining the posterior distribution 
$p(\boldsymbol{r}|\boldsymbol{y}^{(GT)},\boldsymbol{s})$). The highest probability value yields the best estimate for the cylinder position. This process is analogous to a fish attempting to localize the position of a predator or prey. Using the fact that the prior distribution of the disturbance location is uniform, the required posterior probability distribution of the disturbance location is proportional to the likelihood 
$p(\boldsymbol{y}|\boldsymbol{r}^{(GT)},\boldsymbol{s})$ defined in (2.3), where 
$\boldsymbol{y}$ are the measurements recorded along the swimmer body.
The resulting probability distribution for estimating the correct 
$\boldsymbol{r}^{(GT)}$ is depicted in figure 12, with the rows showing results for an increasing number of sensors. The left column shows probability estimates from a self-propelled swimmer attempting to identify the position of an unknown disturbance source, based on flow measurements from optimal sensor locations which are indicated on the body. The right column depicts estimates made by a swimmer using suboptimal sensor distributions. The probability distributions indicate that the swimmer on the left (using optimal sensor distributions) is able to provide a much more accurate estimate of the correct location of the disturbance source. We observe that the uninformed placement of a single sensor in figure 12(b) leads to a large spread in the probability distribution, making it difficult to locate the disturbance source accurately. In comparison, the first optimal sensor in figure 12(a) yields a noticeably narrower spread, centred close to the correct position of the signal-generating cylinder (i.e. the ground-truth). The probability distributions in both cases become narrower with increasing number of sensors, making it easier to locate the disturbance source. In all cases, the optimal arrangement of sensors performs noticeably better than the uniform distribution for identifying the correct cylinder location.
Figure 12. Plots showing the probability that the signal being measured by a self-propelled swimmer originates from a particular location in the prior-region. (a,c,e) Swimmer using flow measurements from optimal sensor locations indicated on the body, panels (b,d,f) show estimates made by a swimmer using suboptimal sensor distributions. Brighter areas indicate regions of higher probability. The relevant sensor arrangement is shown using red ‘
$\times$’ symbols on the swimmers’ bodies. The actual position of the signal-generating cylinder is marked with red diamonds. (a,c,e) Probability distributions computed using measurements from 1, 3 and 5 optimal sensors on both sides of the body. (b,d,f) Probability distributions for 1, 3 and 5 uniformly distributed sensors on both sides.
4 Discussion
While the present work attempts to represent realistic flow conditions, it is important to keep in mind its limitations and the simplifying assumptions used in the simulations and data analyses. Most importantly, we note that flexural dynamics of the sensory structures are not accounted for in the current study. In actual fish, the extent to which the hair-like sensory structures are deflected by the flow can influence their sensing effectiveness (Hudspeth & Corey Reference Hudspeth and Corey1977; Ćurčić-Blake & van Netten Reference Ćurčić-Blake and van Netten2006; Van Trump & McHenry Reference Van Trump and McHenry2008; Bleckmann & Zelick Reference Bleckmann and Zelick2009). For instance, beyond a certain locomotion speed, the superficial neuromasts may suffer from saturation effects during forward motion, when the hair cells are fully deflected and may no longer be able to detect external stimuli effectively. In addition to natural swimmers, such dynamics have been included in artificial lateral lines employed in robotic devices (Fan et al. Reference Fan, Chen, Zou, Bullen, Liu and Delcomyn2002). Similarly, the canal neuromasts are located in recessed channels under the skin, which alter the flow experienced by the sensory structures considerably. This is compensated for in real fish through the introduction of resonant hair cell structures (Ó Maoiléidigh, Nicola & Hudspeth Reference Ó Maoiléidigh, Nicola and Hudspeth2012), but is not accounted for in our simulations. All of these aspects may play an important role in determining the observed distribution of sensors on a fish’s body, in addition to the fluid-flow induced by external stimuli. However, our two-dimensional simulations do not account for these factors, in an attempt to keep the complexity of the Navier–Stokes simulations to manageable levels.
Furthermore, we note that detecting flows by a stationary swimmer in a quiescent fluid with external perturbations is not equivalent to flow measurements by moving swimmers; the latter suffer from separation effects of the boundary layer, which can influence the location of optimal sensors. Consequently, we must be careful when considering sensor placement using the combined datasets, as done in figures 10 and 11. Inspecting the individual scenarios in figures 5(g), 8(b) and 9(b), we observe that the sensor distributions for the stationary and moving swimmers are not entirely dissimilar, which gives us confidence in using the combined dataset for sensor placement in figures 10 and 11.
In our approach, a sensor-pair is placed on the body surface symmetrically around the fish centre line. We anticipate that re-evaluating the simulations with a left–right flip, i.e. changing the orientation of the fish with respect to its approach to the cylinder, may lead to some differences in the exact sensor locations. However, once the fish is experiencing the cylinder’s wake on its pressure and shear sensors, we expect that the overall sensor arrangement will remain unchanged, i.e. we will still observe a dense distribution of sensors in the head and the tail. We have found that the dominant factor in determining sensor placement is not the approach to the flow but rather the body geometry and motion. This is evident when we compare sensor distributions across very dissimilar scenarios, i.e. the three different flow-configurations involving stationary fish, self-propelled fish, and rigid fish placed in the D-cylinder’s wake (figures 5g, 9b and 8b).
The present approach is computationally demanding as it requires conducting a large number of direct numerical simulations (DNS) for the Navier–Stokes equations, each of which takes four to seven hours to complete on a 12-core central processing unit (CPU) node depending on the flow configuration. The computational cost is substantial, given that close to 400 distinct simulations have to be evaluated for each of the considered flow configurations, resulting in a total of approximately 3000 simulations. At the same time, we must emphasize that all DNS computations are performed once (offline) and the sequential sensor placement algorithm is rather inexpensive, with the computations taking of the order of five minutes using a single CPU core.
Hence, it is imperative to store DNS data containing all possible sensor-measurements offline, and to process them for sequential sensor placement as required. We remark that if the same study were to be conducted in a three-dimensional setting using DNS, the computational cost of each simulation would increase by approximately three orders of magnitude, since the computational grid would increase in size from 4096 x 4096 cells for the present two-dimensional (2-D) cases, to approximately 4096 x 4096 x 1024 cells for the three-dimensional (3-D) cases. This would represent a significant increase in computational cost, especially considering the need to run approximately 3000 DNSs to examine all of the cases discussed in the present work. We remark that if, for 3-D flows, we use instead approaches such as large eddy simulations and unsteady Reynolds-averaged Navier–Stokes calculations we do not expect to have higher computational costs than the ones required herein for the 2-D DNS. In closing, we argue that the availability of computational power and automation in the near future will make approaches such as the one presented herein amenable to further computational (and experimental) studies.
5 Conclusion
We have combined 2-D Navier–Stokes simulations with Bayesian optimal experimental design, to identify the best arrangement of sensory structures on self-propelled swimmers’ bodies. The study is inspired by the particular distribution of flow-sensing mechanoreceptors found in many fish species, referred to as the lateral-line organ, where a large number of sensory structures are located in the head.
We optimize sensor arrangements on two different swimmer shapes under the influence of various sources of disturbance. We find optimal arrangements that resemble those found in fish bodies, suggesting that such arrangements may allow them to gather information from their surroundings more effectively than other layouts. We demonstrate that the optimal configuration of these sensors depends on the body shape and the type of disturbance being perceived. This is explored using a variety of simulations involving both static and swimming configurations, using distinct body profiles resembling fish larvae and adults, and using disturbances generated by oscillating cylinders, rotating cylinders and by D-shaped half-cylinders. Despite certain differences that exist in sensor distributions among the various cases considered, there is a marked tendency for a large number of shear stress sensors to be located in the head and the tail of the swimmer, with virtually no sensors found in the midsection. In the case of pressure gradient sensors, we observe a high density of sensors placed in the head, followed by regularly spaced distribution along the entire body. These observations closely reflect the structure of the sensory organ in live fish.
To assess the effectiveness of the sensor placement algorithm, we compare the performance of optimal arrangements to that of uninformed uniform sensor distributions. The results confirm that optimal distribution patterns lead to more accurate identification of external disturbances, which suggests that these distinctive distributions may allow fish to assimilate maximum information from their surroundings using the fewest number of neuromasts. We believe that the present work is a positive step towards understanding mechanosensing in fish, and we hope that the proposed methodology can assist in the development of optimal sensory-layouts for engineered swimmers.
Acknowledgements
This work was supported by European Research Council Advanced Investigator Award 341117, and utilized computational resources granted by the Swiss National Supercomputing Centre (CSCS) under project IDs ‘s658’ and ‘ch7’. We thank P. E. Hadjidoukas for assistance with the TORC interface for running the flow simulations.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2019.940.
Appendix A. Shape parameterisation and swimming kinematics
The adult shape used in the simulations is modelled using three piecewise polynomials
$$\begin{eqnarray}\displaystyle h(s)=\left\{\begin{array}{@{}ll@{}}\displaystyle \sqrt{2sh_{head}-s^{2}},\quad & 0\leqslant s<s_{head},\\ \displaystyle h_{head}-(h_{head}-h_{tail})\left(\frac{s-s_{head}}{s_{tail}-s_{head}}\right)^{2},\quad & s_{head}\leqslant s<s_{tail},\\ \displaystyle h_{parabola},\quad & s_{tail}\leqslant s\leqslant L.\end{array}\right. & & \displaystyle\end{eqnarray}$$ Here, 
$s$ denotes the curvilinear coordinate running from the head to the tail along the body midline, and 
$h$ is the half-width of the body. The constant parameters that determine the final shape are 
$s_{head}=0.04L$, 
$s_{tail}=0.95L$, 
$h_{head}=0.04L$, 
$h_{tail}=0.01L$. The parabolic section that defines the smooth tail is computed as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle a=h_{tail}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle b=\frac{h_{tail}-h_{head}}{s_{tail}-s_{head}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c=\frac{-b(L-s_{tail})-h_{tail}}{(L-s_{tail})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{parabola}=a+b(s-s_{tail})+c(s-s_{tail})^{2}. & \displaystyle\end{eqnarray}$$The larva shape used in the simulations is based on silhouettes of zebrafish extracted from experiments. (The authors thank the Engert Lab at Harvard University for providing the experimental images.) The segmented shape is parameterized using a natural cubic spline comprised of seven piecewise sections, with the following knots and polynomial coefficients:
$$\begin{eqnarray}\displaystyle & \displaystyle (s_{0},\ldots \,s_{i},\ldots ,s_{7})/L=(0.0,0.018,0.058,0.098,0.198,0.238,0.698,1.0),\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a_{i,j}=\left(\begin{array}{@{}cccc@{}}0.000000 & 3.152700 & -44.18100 & 145.49\\ 0.043282 & 1.703600 & -36.32400 & 300.38\\ 0.072530 & 0.239460 & -0.278770 & -54.342\\ 0.078185 & -0.043688 & -6.799800 & 38.155\\ 0.043973 & -0.258990 & 4.646700 & -37.601\\ 0.038641 & -0.067740 & 0.134620 & -0.12819\\ 0.023488 & -0.025271 & -0.042293 & -0.43563\end{array}\right). & \displaystyle\end{eqnarray}$$
$i\in [1,7]$ denotes the polynomial segment between knots 
$s_{i-1}$ and 
$s_{i}$. The corresponding cubic polynomial describing the body half-width in each of the seven sections is given by 
$$\begin{eqnarray}h_{i}(s)=a_{i,1}+a_{i,2}(s-s_{i-1})+a_{i,3}(s-s_{i-1})^{2}+a_{i,4}(s-s_{i-1})^{3}.\end{eqnarray}$$ The swimmers propel themselves by imposing a sinusoidal wave travelling along the body. This wave is specified by the spatially and temporally varying curvature 
$\unicode[STIX]{x1D705}(s,t)$ of the midline,
$$\begin{eqnarray}\unicode[STIX]{x1D705}(s,t)=A(s)\sin \left(\frac{2\unicode[STIX]{x03C0}t}{T}-\frac{2\unicode[STIX]{x03C0}s}{L}+\unicode[STIX]{x1D719}\right),\end{eqnarray}$$ where 
$A(s)$ is described using a natural cubic spline, with six control points located at 
$(s_{i},A_{i})=((0,0.82),(0.15L,1.47),(0.4L,2.57),(0.65L,3.75),(0.9L,5.09),(1L,5.7))$. Here, 
$s$ denotes the curvilinear coordinate running from the head to the tail along the body midline. The time period of body undulation is set to 
$T=0.4$, and the phase difference 
$\unicode[STIX]{x1D719}$ is initially set to 0. The time-varying curvature is imposed as a function of the arclength 
$s$, ensuring that the swimmer’s length remains constant during deformation. The lateral and longitudinal coordinates (as well as the deformation velocities) of the midline points are recovered from the curvature using the Frenet–Serret formulas. The travelling waveform yields steadily swimming fish. However, several species employ an intermittent swimming mode referred to as ‘burst and coast’ swimming (Weihs Reference Weihs1974). The ‘burst and coast’ mode may not only be energetically efficient, but it also allows fish to better accumulate information about external disturbances, by stabilizing the sensory fields (Kramer & McLaughlin Reference Kramer and McLaughlin2001). The ‘burst and coast’ strategy is observed in blind cave fish, where they accelerate and glide past unfamiliar objects and obstacles repeatedly (von Campenhausen et al. Reference von Campenhausen, Riess and Weissert1981). This allows them to form a ‘hydrodynamic image’ of their surroundings, by perceiving reflections of their self-generated motion. Intermittent swimming has also been shown to be critical for avoiding collisions when approaching a wall (Windsor et al. Reference Windsor, Tan and Montgomery2008). The gliding phase during burst-coast swimming may minimize self-generated ‘noise’ in the boundary layer on the body, thereby allowing signals of external origin to permeate through to the neuromasts. Modelling the burst-coast motion for the swimmers involves multiplying the curvature amplitude 
$A(s)$ with a smoothly varying piecewise function 
$f(t)$, as described in Verma et al. (Reference Verma, Hadjidoukas, Wirth and Koumoutsakos2017) as follows:
$$\begin{eqnarray}f(t)=\left\{\begin{array}{@{}ll@{}}1, & t\in \unicode[STIX]{x0394}t_{steady},\\ 1-3\unicode[STIX]{x1D706}_{coast}^{2}+2\unicode[STIX]{x1D706}_{coast}^{3}, & t\in \unicode[STIX]{x0394}t_{decel},\\ 0, & t\in \unicode[STIX]{x0394}t_{coast},\\ 3\unicode[STIX]{x1D706}_{burst}^{2}-2\unicode[STIX]{x1D706}_{burst}^{3}, & t\in \unicode[STIX]{x0394}t_{accel}.\end{array}\right.\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D706}_{coast},\unicode[STIX]{x1D706}_{burst}\in [0,1]$ are ramp functions increasing linearly from 0 to 1 within the intervals 
$\unicode[STIX]{x0394}t_{decel}$ and 
$\unicode[STIX]{x0394}t_{accel}$. The four time intervals describing the burst-coast phases are set to 
$\unicode[STIX]{x0394}t_{steady}=1.5T$, 
$\unicode[STIX]{x0394}t_{decel}=0.375T$, 
$\unicode[STIX]{x0394}t_{coast}=1.0T$ and 
$\unicode[STIX]{x0394}t_{accel}=0.375T$. At the end of each coasting phase, the phase angle in (A 5) is updated to 
$\unicode[STIX]{x1D719}=(\unicode[STIX]{x03C0}-2\unicode[STIX]{x03C0}t_{SA}/T)$, where 
$t_{SA}$ denotes the time at the start of the acceleration phase (or, equivalently, at the end of the coasting phase). This introduces mirror symmetry between tail-beats from one burst-coast cycle to another, and allows the fish to swim on a relatively straight trajectory.
Appendix B. Derivation of the utility estimator
By applying Bayes’ theorem (2.5) can be written equivalently as
$$\begin{eqnarray}U(\boldsymbol{s})=\int _{{\mathcal{Y}}}\int _{{\mathcal{R}}}\ln \frac{p(\boldsymbol{y}|\boldsymbol{r},\boldsymbol{s})}{p(\boldsymbol{y}|\boldsymbol{s})}\;p(\boldsymbol{r})\;p(\boldsymbol{y}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y},\end{eqnarray}$$ where we have used the assumption that 
$p(\boldsymbol{r})=p(\boldsymbol{r}|\boldsymbol{s})$. We approximate the integral over 
${\mathcal{R}}$ with a quadrature rule on the points 
$\boldsymbol{r}^{(i)}$ and weights 
$w_{i}$ for 
$i=1,\ldots ,N_{\boldsymbol{r}}$. Then (B 1) is approximated by
$$\begin{eqnarray}U(\boldsymbol{s})\approx \mathop{\sum }_{i=1}^{N_{\boldsymbol{r}}}w_{i}\;p(\boldsymbol{r}^{(i)})\int _{{\mathcal{Y}}}\ln \frac{p(\boldsymbol{y}|\boldsymbol{r}^{(i)},\boldsymbol{s})}{p(\boldsymbol{y}|\boldsymbol{s})}\;p(\boldsymbol{y}|\boldsymbol{r}^{(i)},\boldsymbol{s})\,\text{d}\boldsymbol{y}.\end{eqnarray}$$ For each 
$i$ in the summation, the integral over 
${\mathcal{Y}}$ is approximated by Monte Carlo integration using 
$N_{\boldsymbol{y}}$ points 
$\{\boldsymbol{y}^{(i,j)}\}_{j=1}^{N_{\boldsymbol{y}}}$ from 
$p_{\boldsymbol{y}}(\,\cdot \,|\boldsymbol{r}^{(i)},\boldsymbol{s})$, leading to
$$\begin{eqnarray}U(\boldsymbol{s})\approx \mathop{\sum }_{i=1}^{N_{\boldsymbol{r}}}\mathop{\sum }_{j=1}^{N_{\boldsymbol{y}}}\frac{w_{i}\;p(\boldsymbol{r}^{(i)})}{N_{\boldsymbol{y}}}\,[\ln p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r}^{(i)},\boldsymbol{s})-\ln p(\boldsymbol{y}^{(i,j)}|\boldsymbol{s})].\end{eqnarray}$$ The quantity 
$p(\boldsymbol{y}^{(i,j)}|\boldsymbol{s})$ is approximated with the same quadrature rule,
$$\begin{eqnarray}\displaystyle p(\boldsymbol{y}^{(i,j)}|\boldsymbol{s}) & = & \displaystyle \int _{{\mathcal{R}}}\,p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r},\boldsymbol{s})\;p(\boldsymbol{r})\,\text{d}\boldsymbol{r}\nonumber\\ \displaystyle & {\approx} & \displaystyle \mathop{\sum }_{k=1}^{N_{\boldsymbol{r}}}w_{k}\;p(\boldsymbol{y}^{(i,j)}|\boldsymbol{r}^{(k)},\boldsymbol{s})\;p(\boldsymbol{r}^{(k)}).\end{eqnarray}$$Finally, substituting (B 4) into (B 3), the estimator for 
$U(\boldsymbol{s})$ is given by (2.6).
Appendix C. Utility of combined experiments
We assume that the measurements 
$\boldsymbol{y}_{1}$ and 
$\boldsymbol{y}_{2}$ are independent conditioned on 
$\boldsymbol{s}$. We want to show that the utility function using both experiments 
$U_{1,2}$ is equal to the sum of the individual utility functions, 
$U_{1}$ and 
$U_{2}$.
$$\begin{eqnarray}\displaystyle U_{1,2}(\boldsymbol{s}) & = & \displaystyle \int _{{\mathcal{Y}}_{1}}\int _{{\mathcal{Y}}_{2}}\int _{\unicode[STIX]{x1D6E9}}\ln \frac{p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})}{p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{s})}\;p(\boldsymbol{r})\;p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{1}\,\text{d}\boldsymbol{y}_{2}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\mathcal{Y}}_{1}}\int _{{\mathcal{Y}}_{2}}\int _{\unicode[STIX]{x1D6E9}}\ln \frac{p(\boldsymbol{y}_{1}|\boldsymbol{r},\boldsymbol{s})p(\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})}{p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{s})}\;p(\boldsymbol{r})\;p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{1}\,\text{d}\boldsymbol{y}_{2}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\mathcal{Y}}_{1}}\int _{\unicode[STIX]{x1D6E9}}\ln p(\boldsymbol{y}_{1}|\boldsymbol{r},\boldsymbol{s})\,p(\boldsymbol{r})~p(\boldsymbol{y}_{1}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{1}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{{\mathcal{Y}}_{2}}\int _{\unicode[STIX]{x1D6E9}}\ln p(\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})\,p(\boldsymbol{r})~p(\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{2}-h(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{s}),\end{eqnarray}$$where
$$\begin{eqnarray}h(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{s})=\int _{{\mathcal{Y}}_{1}}\int _{{\mathcal{Y}}_{2}}\int _{\unicode[STIX]{x1D6E9}}\ln p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{s})\;p(\boldsymbol{r})\;p(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{1}\,\text{d}\boldsymbol{y}_{2}.\end{eqnarray}$$Using (B 1), it is easy to check that
$$\begin{eqnarray}U_{i}(\boldsymbol{s})=\int _{{\mathcal{Y}}_{i}}\int _{\unicode[STIX]{x1D6E9}}\ln p(\boldsymbol{y}_{i}|\boldsymbol{r},\boldsymbol{s})\,p(\boldsymbol{r})~p(\boldsymbol{y}_{i}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{i}-h(\boldsymbol{y}_{i}|\boldsymbol{s}),\end{eqnarray}$$ for 
$i=1,2$, where 
$h(\boldsymbol{y}_{i}|\boldsymbol{s})$ is given by
$$\begin{eqnarray}h(\boldsymbol{y}_{i}|\boldsymbol{s})=\int _{{\mathcal{Y}}}\int _{{\mathcal{R}}}\ln p(\boldsymbol{y}_{i}|\boldsymbol{s})\;p(\boldsymbol{r})\;p(\boldsymbol{y}_{i}|\boldsymbol{r},\boldsymbol{s})\,\text{d}\boldsymbol{r}\,\text{d}\boldsymbol{y}_{i}.\end{eqnarray}$$Substituting (C 3) into (C 1),
$$\begin{eqnarray}U_{1,2}(\boldsymbol{s})=U_{1}(\boldsymbol{s})+U_{2}(\boldsymbol{s})-[h(\boldsymbol{y}_{1},\boldsymbol{y}_{2}|\boldsymbol{s})-h(\boldsymbol{y}_{1}|\boldsymbol{s})-h(\boldsymbol{y}_{2}|\boldsymbol{s})].\end{eqnarray}$$ Using (C 2) and (C 4), along with the fact that 
$\boldsymbol{y}_{1}$ is independent of 
$\boldsymbol{y}_{2}$ conditioned on 
$\boldsymbol{s}$, the term inside the brackets is equal to zero, leading to
$$\begin{eqnarray}U_{1,2}(\boldsymbol{s})=U_{1}(\boldsymbol{s})+U_{2}(\boldsymbol{s}).\end{eqnarray}$$ The assumption of conditional independency is valid in our case since knowledge of the measurements of one experiment 
$\boldsymbol{y}_{1}$ and the location of the sensors 
$\boldsymbol{s}$ does not provide any additional information for the measurements of the other experiment compared to the case of knowledge only of the sensor location. Equivalently, 
$p(\boldsymbol{y}_{2}|\boldsymbol{y}_{1},\boldsymbol{s})=p(\boldsymbol{y}_{2}|\boldsymbol{s})$.
