Optimal sensor placement for artificial swimmers

TL;DR

This paper investigates the optimal placement of sensors on artificial swimmers to maximize environmental information gathering, revealing arrangements similar to natural fish neuromasts and suggesting evolutionary optimality.

Contribution

It introduces a method combining Navier-Stokes simulations and Bayesian design to identify optimal sensor configurations for flow disturbance detection.

Findings

Shear sensors are best placed at the head and tail.

Pressure sensors are densely located at the head and evenly along the body.

Optimal sensor arrangements resemble natural fish neuromast distributions.

Abstract

Natural swimmers rely for their survival on sensors that gather information from the environment and guide their actions. The spatial organization of these sensors, such as the visual fish system and lateral line, suggests evolutionary selection, but their optimality remains an open question. Here, we identify sensor configurations that enable swimmers to maximize the information gathered from their surrounding flow field. We examine two-dimensional, self-propelled and stationary swimmers that are exposed to disturbances generated by oscillating, rotating and D-shaped cylinders. We combine simulations of the Navier-Stokes equations with Bayesian experimental design to determine the optimal arrangements of shear and pressure sensors that best identify the locations of the disturbance-generating sources. We find a marked tendency for shear stress sensors to be located in the head and the tail of the swimmer, while they are absent from the midsection. In turn, we find a high density of pressure sensors in the head along with a uniform distribution along the entire body. The resulting optimal sensor arrangements resemble neuromast distributions observed in fish and provide evidence for optimality in sensor distribution for natural swimmers.

Tables

Table 1. Table 1 : Optimal distribution of the first 10 shear stress sensors for the self-propelled swimmers. The body has been divided into 3 distinct segments: the head ( 0 ≤ s / L < 0.2 0 𝑠 𝐿 0.2 0\leq s/L<0.2 ); the midsection ( 0.2 ≤ s / L < 0.6 0.2 𝑠 𝐿 0.6 0.2\leq s/L<0.6 ); and the posterior ( 0.6 ≤ s / L ≤ 1 0.6 𝑠 𝐿 1 0.6\leq s/L\leq 1 ).

		Head	Midsection	Posterior
Larva	Oscillating	$1,2,4,5,8,9$	—	$3,6,7,10$
	Rotating	$1,2,\dots ,10$	—	—
Adult	Oscillating	$1,2,4,6,9,10$	—	$3,5,7,8$
	Rotating	$1,2,\dots ,10$	—	—

Equations

$$\frac{\partial \omega }{\partial t}+(\mathbf{u}\cdot \nabla )\omega =\nu {\nabla }^2\omega +\lambda \nabla \times \left(\chi \left({\mathbf{u}}_s-\mathbf{u}\right)\right),$$

$$\mathbf{y}=\mathbf{F}(\mathbf{r};\mathbf{s})+\mathbf{\epsilon }(\mathbf{s}).$$

$$p\left(\mathbf{y}|\mathbf{r},\mathbf{s}\right)$$

$$u(\mathbf{s},\mathbf{y}):={\int }_{\mathcal{R}}p(\mathbf{r}|\mathbf{y},\mathbf{s})\ln \frac{p(\mathbf{r}|\mathbf{y},\mathbf{s})}{p(\mathbf{r})}\text{ d}\mathbf{r}.$$

$$\begin{split}U(\bm{s}):=\mathbb{E}_{\bm{y}|\bm{s}}\big{[}u(\bm{s},\bm{y})\big{]}&=\int_{\mathcal{Y}}u(\bm{s},\bm{y})\;p(\bm{y}|\bm{s}){\text{ d}{\bm{y}}}\\ &=\int_{\mathcal{Y}}\int_{\mathcal{R}}p({\bm{r}}|\bm{y},\bm{s})\;\ln\frac{p({\bm{r}}|\bm{y},\bm{s})}{p({\bm{r}})}\;p(\bm{y}|\bm{s}){\text{ d}{{\bm{r}}}}{\text{ d}{\bm{y}}}\,.\end{split}$$

$$\hat{U}(\mathbf{s})=\frac{1}{N_{\mathbf{y}}}\sum _{j=1}^{N_{\mathbf{y}}}\sum _{i=1}^{N_{\mathbf{r}}}{w}_{i}p({\mathbf{r}}^{(i)})\left[\ln p({\mathbf{y}}^{(i,j)}|{\mathbf{r}}^{(i)},\mathbf{s})-\ln \left(\sum _{k=1}^{N_{\mathbf{r}}}{w}_{k}p({\mathbf{r}}^{(k)})p({\mathbf{y}}^{(i,j)}|{\mathbf{r}}^{(k)},\mathbf{s})\right)\right].$$

$${\mathbf{y}}^{(i,j)}=\mathbf{F}({\mathbf{r}}^{(i)};\mathbf{s})+{\mathbf{\epsilon }}^{(j)},$$

$${\Sigma }_{ij}(\mathbf{s})=\{\begin{array}{ll}{\sigma }^2\exp \left(-\frac{\Vert \mathbf{x}(s_i)-\mathbf{x}(s_j)\Vert }{\ell }\right),& \text{ if }1\le i,j\le n,\\ {\Sigma }_{i-n,j-n}(\mathbf{s}),& \text{ if }n<i,j\le 2n,\\ 0& \text{ otherwise,}\end{array}$$

$$s_i^{\star }=\underset{s}{\mathrm{argmax}}{\hat{U}}_i(s)\text{and}{\hat{U}}_i^{\star }=\underset{s}{\max }{\hat{U}}_i(s).$$

$$h(s)=\{\begin{array}{ll}\sqrt{2s{h}_{head}-s^2}& ,0\le s<s_{head}\\ h_{head}-(h_{head}-h_{tail}){\left(\frac{s-s_{head}}{s_{tail}-s_{head}}\right)}^2& ,s_{head}\le s<s_{tail}\\ h_{parabola}& ,s_{tail}\le s\le L.\end{array}$$

$$a$$

$$b$$

$$c$$

$$h_{parabola}$$

$$\left(s_0,\cdots s_i,\cdots ,s_7\right)/L$$

$$a_{i,j}$$

$$h_i(s)=a_{i,1}+a_{i,2}\left(s-s_{i-1}\right)+a_{i,3}{\left(s-s_{i-1}\right)}^2+a_{i,4}{\left(s-s_{i-1}\right)}^3.$$

$$\kappa (s,t)=A(s)\sin \left(\frac{2\pi t}{T}-\frac{2\pi s}{L}+\varphi \right),$$

$$f(t)=\left\{\begin{array}[]{ll}1&t\in\Delta t_{steady}\\ 1-3\lambda_{coast}^{2}+2\lambda_{coast}^{3}&t\in\Delta t_{decel}\\ 0&t\in\Delta t_{coast}\\ 3\lambda_{burst}^{2}-2\lambda_{burst}^{3}&t\in\Delta t_{accel}\,.\end{array}\right.$$

$$U(\mathbf{s})={\int }_{\mathcal{Y}}{\int }_{\mathcal{R}}\ln \frac{p(\mathbf{y}|\mathbf{r},\mathbf{s})}{p(\mathbf{y}|\mathbf{s})}p(\mathbf{r})p(\mathbf{y}|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}\mathbf{y},$$

$$U(\mathbf{s})\approx \sum _{i=1}^{N_{\mathbf{r}}}{w}_{i}p({\mathbf{r}}^{(i)}){\int }_{\mathcal{Y}}\ln \frac{p(\mathbf{y}|{\mathbf{r}}^{(i)},\mathbf{s})}{p(\mathbf{y}|\mathbf{s})}p(\mathbf{y}|{\mathbf{r}}^{(i)},\mathbf{s})\text{ d}\mathbf{y}.$$

$$U(\bm{s})\approx\sum_{i=1}^{N_{{\bm{r}}}}\sum_{j=1}^{N_{\bm{y}}}\frac{w_{i}\;p({\bm{r}}^{(i)})}{N_{\bm{y}}}\,\Big{[}\ln{p(\bm{y}^{(i,j)}|{\bm{r}}^{(i)},\bm{s})}-\ln{p(\bm{y}^{(i,j)}|\bm{s})}\Big{]}\,.$$

$$\begin{array}{rl}p({\mathbf{y}}^{(i,j)}|\mathbf{s})& ={\int }_{\mathcal{R}}p({\mathbf{y}}^{(i,j)}|\mathbf{r},\mathbf{s})p(\mathbf{r})\text{ d}\mathbf{r}\\ & \approx \sum _{k=1}^{N_{\mathbf{r}}}{w}_{k}p({\mathbf{y}}^{(i,j)}|{\mathbf{r}}^{(k)},\mathbf{s})p({\mathbf{r}}^{(k)}).\end{array}$$

$$\begin{array}{rl}{U}_{1,2}(\mathbf{s})& ={\int }_{{\mathcal{Y}}_1}{\int }_{{\mathcal{Y}}_2}{\int }_{\Theta }\ln \frac{p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{r},\mathbf{s})}{p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{s})}p(\mathbf{r})p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_1\text{ d}{\mathbf{y}}_2\\ & ={\int }_{{\mathcal{Y}}_1}{\int }_{{\mathcal{Y}}_2}{\int }_{\Theta }\ln \frac{p({\mathbf{y}}_1|\mathbf{r},\mathbf{s})p({\mathbf{y}}_2|\mathbf{r},\mathbf{s})}{p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{s})}p(\mathbf{r})p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_1\text{ d}{\mathbf{y}}_2\\ & ={\int }_{{\mathcal{Y}}_1}{\int }_{\Theta }\ln p({\mathbf{y}}_1|\mathbf{r},\mathbf{s})p(\mathbf{r})p({\mathbf{y}}_1|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_1\\ & +{\int }_{{\mathcal{Y}}_2}{\int }_{\Theta }\ln p({\mathbf{y}}_2|\mathbf{r},\mathbf{s})p(\mathbf{r})p({\mathbf{y}}_2|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_2-h({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{s}),\end{array}$$

$$h({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{s})={\int }_{{\mathcal{Y}}_1}{\int }_{{\mathcal{Y}}_2}{\int }_{\Theta }\ln p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{s})p(\mathbf{r})p({\mathbf{y}}_1,{\mathbf{y}}_2|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_1\text{ d}{\mathbf{y}}_2,$$

$$U_i(\mathbf{s})={\int }_{{\mathcal{Y}}_i}{\int }_{\Theta }\ln p({\mathbf{y}}_i|\mathbf{r},\mathbf{s})p(\mathbf{r})p({\mathbf{y}}_i|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_i-h({\mathbf{y}}_i|\mathbf{s}),$$

$$h({\mathbf{y}}_i|\mathbf{s})={\int }_{\mathcal{Y}}{\int }_{\mathcal{R}}\ln p({\mathbf{y}}_i|\mathbf{s})p(\mathbf{r})p({\mathbf{y}}_i|\mathbf{r},\mathbf{s})\text{ d}\mathbf{r}\text{ d}{\mathbf{y}}_i,$$

$$U_{1,2}(\bm{s})=U_{1}(\bm{s})+U_{2}(\bm{s})-\Big{[}h(\bm{y}_{1},\bm{y}_{2}|\bm{s})-h(\bm{y}_{1}|\bm{s})-h(\bm{y}_{2}|\bm{s})\Big{]}\,.$$

$$U_{1,2}(\mathbf{s})=U_1(\mathbf{s})+U_2(\mathbf{s}).$$

Full text

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, 1975; Ladich & Bass, 2003; 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, 1974) 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, 1963; Kroese & Schellart, 1992; Coombs et al., 1996; Coombs & Netten, 2005; Bleckmann & Zelick, 2009). 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. (2000)) or non-zero pressure gradients (sub-surface ‘canal’ neuromasts - Bleckmann & Zelick (2009)). An array of such sensors allows fish to discern both the direction and speed of disturbances generated in their surrounding flow (Chambers et al., 2014; Asadnia et al., 2015).

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 et al., 2015), 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 (Engelmann et al., 2000; Bleckmann & Zelick, 2009; Atema et al., 1988), 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 et al., 1976; Partridge & Pitcher, 1980), predator evasion (Blaxter & Fuiman, 1989), prey detection/capture (Hoekstra & Janssen, 1985), reproduction (Satou et al., 1994), rheotaxis (Dijkgraaf, 1963; Kanter & Coombs, 2003), obstacle avoidance (Hassan, 1989), and station-keeping by countering the effects of unsteady gusts (Sutterlin & Waddy, 1975). Disrupting the normal functioning of the lateral line, either via chemical or mechanical means, hinders fish’s ability to perform these tasks effectively. Liao (2006) demonstrated that disabling the lateral line system influences 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, 1989; Coombs & Montgomery, 1999; Montgomery et al., 2001; Coombs & Braun, 2003). 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 et al., 1981; Windsor et al., 2008; de Perera, 2004).

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., 2006, 2010; Kottapalli et al., 2012; Ježov et al., 2012; Kruusmaa et al., 2014; Asadnia et al., 2015; Triantafyllou et al., 2016; Strokina et al., 2016; Kottapalli et al., 2018; Yen et al., 2018). 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., 2011; Venturelli et al., 2012; Xu & Mohseni, 2017), and mechanically deforming hair-like structures (Yang et al., 2006; Tao & Yu, 2012; Abdulsadda & Tan, 2013; Dagamseh et al., 2013; DeVries et al., 2015; Triantafyllou et al., 2016).

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, 1984; Denton & Gray, 1988; Kroese & Schellart, 1992; Coombs et al., 1996; Ćurčić-Blake & van Netten, 2006). 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, 1998), and Karman vortex streets (Venturelli et al., 2012). 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, 1992; Franosch et al., 2009; Bouffanais et al., 2011; Ren & Mohseni, 2012; Colvert & Kanso, 2016). A few of these studies have attempted to infer the optimal arrangement of sensors on rigid objects exposed to various flow conditions. Colvert & Kanso (2016) determined the optimal placement of a single sensor-pair on an elliptical body, moving at different orientations in uniform flow. Ahrari et al. (2017) 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., 2016), 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, 2014), which cannot be accounted-for via inviscid assumptions. The importance of viscous effects has also been demonstrated by Rapo et al. (2009), 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 (Colvert et al., 2018; Alsalman et al., 2018). 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,

[TABLE]

where $\bm{u}$ is the flow-velocity, and $\omega=\nabla\times{\bm{u}}$ is the vorticity. The penalty term, $\lambda\nabla\times\left(\chi\left(\bm{u}_{s}-\bm{u}\right)\right)$ models the interaction of objects with the surrounding fluid (Coquerelle & Cottet (2008)), where $0<\chi\leq 1$ indicates the solid body. $\lambda$ is the penalization parameter, and $\bm{u}_{s}$ represents the combined translational, rotational, and deformational velocity of the solid object. The equations are discretized using remeshed vortex methods (Koumoutsakos & Leonard, 1995) and wavelet adapted grids (Rossinelli et al., 2015), 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. (2011) and Rossinelli et al. (2015). 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 mid-line.

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 2(b)) (Gazzola et al., 2012). Details regarding shape parametrization for both cases are provided in the Appendix (Eqs. 10 and 12).

The swimmers propel themselves by imposing a sinusoidal wave travelling along the body. Details of the swimming kinematics are also provided in Appendix A.

2.3 Disturbance-generation and detection

The sensory cues detected by the rigid and swimming bodies described in section 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}=10Hz$, whereas the angular velocity of the rotating cylinders is set to $20\pi\ rad/s$ (10 rotations/s). The cylinders are placed at various locations within a prescribed region in the computational domain (in the ‘prior-region’), as shown in Figure 3.

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 2 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, 1992; Bleckmann & Zelick, 2009; Asadnia et al., 2015). 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 [math]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 section 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 ${\bm{r}}=\left(x,y\right)$ 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 (${\bm{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 $\bm{s}\in R^{n}$, that is the mid-line 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 $\bm{s}$, and are listed in a vector $\bm{y}\in R^{2n}$.

We denote as $\bm{F}({\bm{r}};\bm{s})$ the predictions of shear stress /pressure gradient at sensor locations $\bm{s}$, obtained by solving the Navier-Stokes equations with a disturbance-generating source located at ${\bm{r}}$. Moreover, we assume that we have prior knowledge about the parameter ${\bm{r}}$, encoded in a prior probability distribution $p({\bm{r}})$.

After observing the measurements $\bm{y}$ from sensors $\bm{s}$, we use Bayesian inference to update our prior belief for the plausible values of parameter ${\bm{r}}$, by identifying the posterior probability distribution $p({\bm{r}}|\bm{y},\bm{s})$. Following Bayes’ rule, the posterior distribution $p({\bm{r}}|\bm{y},\bm{s})$ of the model parameters is proportional to the product of the prior distribution $p({\bm{r}})$ and the likelihood $p(\bm{y}|{\bm{r}},\bm{s})$. The likelihood function represents the probability that a particular measurement $\bm{y}$ for a given sensor arrangement $\bm{s}$ originates from the disturbance-source located at ${\bm{r}}$. We assume a prediction error, $\bm{\varepsilon}(\bm{s})$, as the difference between the measurements $\bm{y}$ and the predictions $\bm{F}({\bm{r}};\bm{s})$ such that:

[TABLE]

The prediction error term ($\bm{\varepsilon}(\bm{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 $\bm{\varepsilon}(\bm{s})$ follows a multivariate Gaussian distribution $\mathcal{N}(0,\Sigma(\bm{s}))$ with zero mean and covariance matrix $\Sigma(\bm{s})\in R^{2n\times 2n}$. The likelihood function $p(\bm{y}|{\bm{r}},\bm{s})$ is then expressed as:

[TABLE]

2.4.2 Optimal sensor placement based on information gain

The goal of the optimal sensor placement problem is to find the locations $\bm{s}$ of the sensors such that the data measured in these locations are most informative for estimating the position ${\bm{r}}$ of the disturbance. A measure of information gain is provided by the Kullback-Leibler (KL) 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:

[TABLE]

We note that in the experimental design phase, the measurements $\bm{y}$ are not available. Thus, the prediction error model (Eq. 2) is used to generate measurements for given model parameter values ${\bm{r}}$ and sensor configuration $\bm{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 Eq. 2 (Ryan, 2003):

[TABLE]

The expected utility function involves a double integral over the parameter space ${\bm{r}}$ and over the measured data $\bm{y}$. An efficient estimator of this double integral using sampling techniques is provided by Huan & Marzouk (2013). A similar estimator is used in the present work,

[TABLE]

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 ${\bm{r}}$. In Eq. 6, ${\bm{r}}^{(i)}$ and $w_{i}$ denote $N_{{\bm{r}}}$ quadrature points and corresponding weights related to discretization of the two-dimensional prior-region (Figure 3). A total of $N_{{\bm{r}}}$ distinct Navier-Stokes simulations are conducted, with a cylinder positioned at various discrete points ${\bm{r}}^{(i)}$, and the quadrature is evaluated using the trapezoidal rule. Based on the prediction error defined in Eq. 2, the measured data $\bm{y}^{(i,j)}$ in Eq. 6 are given by,

[TABLE]

where \bm{\varepsilon}^{(j)},\,{\color[rgb]{0,0,0}\text{ with }}j=1,\ldots,N_{\bm{y}}, are vectors sampled from the distribution $\mathcal{N}(0,\Sigma(\bm{s}))$. $N_{\bm{y}}$ is set to $100$ in the current work, which results in a smoother estimate of $\hat{U}(\bm{s})$ in Eq. 6.

We note that the computational effort for evaluating $\hat{U}(\bm{s})$ in Eq. 6 depends primarily on the number of Navier-Stokes simulations, $N_{{\bm{r}}}$, which are required to evaluate $\bm{F}({\bm{r}}^{(i)},\bm{s})$ for different disturbance locations ${\bm{r}}^{(i)}$, and subsequently to determine $\bm{y}^{(i,j)}$ using Eq. 7. The computational burden does not depend on the number of measured samples $N_{\bm{y}}$, since there are no additional time consuming simulations involved in generating $\bm{\varepsilon}^{(j)}$. Thus the computational effort scales linearly with the number $N_{r}$ of model parameter points ${\bm{r}}^{(i)}$.

We assume that the prior distribution $p({\bm{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 ${\bm{r}}$, $p({\bm{r}}|\bm{y},\bm{s})$, is proportional to the likelihood function $p(\bm{y}|{\bm{r}},\bm{s})$.

The covariance matrix $\Sigma(\bm{s})$ depends primarily on the sensor positions $\bm{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,

[TABLE]

where $\bm{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 $\sigma$ 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 $\sigma$ is a fixed percentage ($30\%$) of the mean sensor-measurement, which is computed over all available instances of ${\bm{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, 2012; Simoen et al., 2013), and prevents excessive clustering of sensors within confined neighbourhoods.

Finally, we provide an intuitive interpretation of how Eq. 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 $\bm{y}^{(i,j)}$ has been generated by a disturbance located at ${\bm{r}}^{(i)}$. This implies that the posterior $p({\bm{r}}|\bm{y}^{(i,j)},\bm{s})$, which indicates the probability that a particular disturbance source ${\bm{r}}$ has generated the measurement measurement $\bm{y}^{(i,j)}$, is peaked and centered around the true source location ${\bm{r}}^{(i)}$. Since the prior distribution is uniform, the likelihood $p(\bm{y}^{(i,j)}|{\bm{r}},\bm{s})$ is proportional to the posterior, and is also peaked and centered around ${\bm{r}}^{(i)}$. Thus, the first term in Eq. 6 is large, whereas most of the terms in the second sum are close to zero, since the probability of measurement $\bm{y}^{(i,j)}$ originating from source ${\bm{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 Eq. 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 Eq. 6 larger (non-zero $p(\bm{y}^{(i,j)}|{\bm{r}}^{(k)},\bm{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}(\bm{s})$ described in Eq. 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 (2004), have been demonstrated to be effective alternatives. In this approach, the optimization is carried out iteratively, one sensor at a time. First, $\hat{U}(\bm{s})$ is computed for a single sensor-pair $\bm{s}=s_{1}$, and the optimal solution $s_{1}^{\star}$ is obtained by identifying the maximum in $\hat{U}(\bm{s})$. Then, $\hat{U}(\bm{s})$ is recomputed with $\bm{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,

[TABLE]

We note that the scalar variable $s$ denotes the position of a single sensor-pair, whereas the vector $\bm{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 $\bm{s}^{\star}=(s_{1}^{\star},\ldots,s_{N_{s}}^{\star})$. Papadimitriou (2004) 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\Delta s,\;k=0,\ldots,N_{g}\,\}$, with $\Delta s=L/N_{g}$ and $N_{g}=1000$. Thus, for each iteration of sequential optimization, the utility estimator in Eq. 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 $\bm{s}$. To minimize the relevant computational cost, we run $N_{{\bm{r}}}$ distinct Navier-Stokes simulations for all disturbance locations ${\bm{r}}^{(i)}$ ($i=1\ldots,N_{r}$), and store the shear stress and pressure gradient at all available dicretization 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$.

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 setup 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 setup is shown in Figure 4, and depicts the vorticity generated by the cylinder along the larva’s body.

These two setups allow us to analyze 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 setups corresponds to ${\bm{r}}_{min}=(0.357,0.375)$ and ${\bm{r}}_{max}=(0.7,0.47)$. A total of $11\times 37=407$ potential ${\bm{r}}^{(i)}$ locations are distributed uniformly throughout the region, and the static object’s center of mass is located at $(0.5,0.3)$. The kinematic viscosity is set to $\nu=1e\text{-}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}(\bm{s})$ defined in Eq. 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) to 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 5(a)). The utility curve exhibits consistently high values for $s/L\leq 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 5(e)), with almost no variation along the midsection ($0.2<s/L\leq 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\leq s/L\leq 0.6$ in Figure 5(c)), 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 9 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\leq 0.15$, Figure 5(b)) correspond to large deviations in both $u$ and $v$ near the head (Figures 5(d) and 5(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.

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 section 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 5(g)), 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 Eq. 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.  (Liao et al., 2003) who examined the fluid dynamics of trout placing themselves behind rocks. A uniform horizontal flow of $1L/s$ 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 ${\bm{r}}_{min}=(0.3,0.43)$ to ${\bm{r}}_{max}=(0.79,0.57)$. A total of $11\times 36=396$ potential ${\bm{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 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\leq s/L\leq 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 section 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.

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, 2001), 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 gain a speed of approximately {\color[rgb]{0,0,0}0.7}L/s, which corresponds to a Reynolds number of \text{Re}=uL/\nu\approx{\color[rgb]{0,0,0}280} (with $L=0.2$ and $\nu=1e\text{-}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 ${\bm{r}}_{min}=(0.25,0.375)$ to ${\bm{r}}_{max}=(0.7,0.5)$, with a total of $11\times 37=407$ potential ${\bm{r}}^{(i)}$ locations distributed uniformly throughout the region. The swimmer’s center 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\leq 0.2$, which corresponds to large curvature changes in the body surface. The utility curve also shows a gradual variation for $s/L\geq 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\leq 0.05$ and more gradually for $s/L\geq 0.6$. Furthermore, we note that the blue utility curves in Figures 9(a) and 9(c) are close to [math] for $s/L\geq 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. 9(a)  9(c)

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 1(a)), where large numbers are found in the head and the tail, with sparser clustering in the midsection.

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 setups 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.

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., 2015). The canal is usually present in a continous line running from head to tail, and shows a high concentration of neuromasts in canal branches found in the head (Coombs et al., 1988; Ristroph et al., 2015). 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 section 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 11(b)), and that of shear stress sensors (Figure 10(b)), 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\leq 0.1$, which corresponds to high utility values in Figure 11(a)), 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. (2015), 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. (2015) 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 $\bm{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 $\bm{y}^{(GT)}$ at the sensor locations, generated from a single cylinder located at ${\bm{r}}^{(GT)}$ (the superscript $GT$ denotes ‘ground-truth’). For a given sensor configuration $\bm{s}$, the measurements $\bm{y}^{(GT)}$ are computed using the prediction error model $\bm{y}^{(GT)}=\bm{F}({\bm{r}}^{(GT)};\bm{s})+\bm{\varepsilon}(\bm{s})$, where $\bm{F}({\bm{r}}^{(GT)};\bm{s})$ is obtained by simulating the Navier-Stokes equations with an oscillating cylinder located at ${\bm{r}}^{(GT)}$, and $\bm{\varepsilon}(\bm{s})$ is a vector sampled from the Gaussian distribution $\mathcal{N}(0,\Sigma(\bm{s}))$. Assuming that the disturbance position ${\bm{r}}^{(GT)}$ is unknown, the swimmer attempts to identify it by assigning probability values to all possible cylinder locations ${\bm{r}}$ within the prior-region (i.e., by determining the posterior distribution $p({\bm{r}}|\bm{y}^{(GT)},\bm{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(\bm{y}|{\bm{r}}^{(GT)},\bm{s})$ defined in Eq. 3, where $\bm{y}$ are the measurements recorded along the swimmer body.

The resulting probability distribution for estimating the correct ${\bm{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 un-informed 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, centered 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.

4 Discussion

While the present work attempts to represent realistic flow conditions, it is important to keep in mind its limitations and 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, 1977; Ćurčić-Blake & van Netten, 2006; Van Trump & McHenry, 2008; Bleckmann & Zelick, 2009). 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., 2002). 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 in real fish through the introduction of resonant hair cell structures (Ó Maoiléidigh et al., 2012), 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 centerline. 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 5(g), 9(b), and 8(b)).

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 4 to 7 hours to complete on a 12-core 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 (off-line) and the sequential sensor-placement algorithm is rather inexpensive, with the computations taking on the order of 5 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 Direct Numerical Simulations, 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 2D cases, to approximately 4096 x 4096 x 1024 cells for the 3D 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 3D flows, we use instead approaches such as Large Eddy Simulations (LES) and Unsteady Reynolds Averaged Navier-Stokes (URANS) calculations we do not expect to have higher computational costs than the ones required herein for the 2D 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 two-dimensional 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 un-informed 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.

6 Acknowledgments

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 Panagiotis E. Hadjidoukas for assistance with the TORC interface for running the flow simulations.

Appendix A Shape parametrisation and swimming kinematics

The adult shape used in the simulations is modelled using three piecewise polynomials,

[TABLE]

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:

[TABLE]

The larva shape used in the simulations is based on silhouettes of zebrafish extracted from experiments111The authors thank the Engert Lab at Harvard University for providing the experimental images.. The segmented shape is parametrized using a natural cubic spline comprised of 7 piecewise sections, with the following knots and polynomial coefficients:

[TABLE]

Here, subscript $i\in\left[1,7\right]$ 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 7 sections is given by

[TABLE]

The swimmers propel themselves by imposing a sinusoidal wave travelling along the body. This wave is specified by the spatially- and temporally-varying curvature $\kappa(s,t)$ of the mid-line,

[TABLE]

where $A(s)$ is described using a natural cubic spline, with 6 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 mid-line. The time period of body-undulation is set to $T=0.4$, and the phase difference $\phi$ is initially set to 0. The time-varying curvature is imposed as a function of the arc-length $s$, ensuring that the swimmer’s length remains constant during deformation. The lateral and longitudinal coordinates (as well as the deformation velocities) of the mid-line points are recovered from the curvature using the Frenet-Serret formulas. The travelling wave form yields steadily-swimming fish. However, several species employ an intermittent swimming mode referred to as ‘burst-and-coast’ swimming (Weihs, 1974). 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, 2001). 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., 1981). 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., 2008). 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. Modeling the 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. (2017):

[TABLE]

Here, $\lambda_{coast},\lambda_{burst}\in[0,1]$ are ramp functions increasing linearly from 0 to 1 within the intervals $\Delta t_{decel}$ and $\Delta t_{accel}$. The 4 time-intervals describing the burst-coast phases are set to $\Delta t_{steady}=1.5T$, $\Delta t_{decel}=0.375T$, $\Delta t_{coast}=1.0T$, and $\Delta t_{accel}=0.375T$. At the end of each coasting phase, the phase angle in Eq. 14 is updated to $\phi=\left(\pi-2\pi t_{SA}/T\right)$, 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 Eq. 5 can be written equivalently as,

[TABLE]

where we have used the assumption that $p({\bm{r}})=p({\bm{r}}|\bm{s})$. We approximate the integral over $\mathcal{R}$ with a quadrature rule on the points ${\bm{r}}^{(i)}$ and weights $w_{i}$ for $i=1,\ldots,N_{{\bm{r}}}$. Then Eq. 16 is approximated by,

[TABLE]

For each $i$ in the summation, the integral over $\mathcal{Y}$ is approximated by Monte Carlo integration using $N_{\bm{y}}$ points $\{\bm{y}^{(i,j)}\}_{j=1}^{N_{\bm{y}}}$ from $p_{\bm{y}}(\,\cdot\,|{\bm{r}}^{(i)},\bm{s})$, leading to

[TABLE]

The quantity $p(\bm{y}^{(i,j)}|\bm{s})$ is approximated with the same quadrature rule,

[TABLE]

Finally, substituting Eq. 19 into Eq. 18, the estimator for $U(\bm{s})$ is given by Eq. 6.

Appendix C Utility of combined experiments

We assume that the measurements $\bm{y}_{1}$ and $\bm{y}_{2}$ are independent conditioned on $\bm{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}$.

[TABLE]

where

[TABLE]

Using Eq. 16, it is easy to check that

[TABLE]

for $i=1,2$, where $h(\bm{y}_{i}|\bm{s})$ is given by

[TABLE]

Substituting Eq. 22 into Eq. 20,

[TABLE]

Using Eq. 21 and Eq. 23, along with the fact that $\bm{y}_{1}$ is independent of $\bm{y}_{2}$ conditioned on $\bm{s}$, the term inside the brackets is equal to zero, leading to

[TABLE]

The assumption of conditional independency is valid in our case since knowledge of the measurements of one experiment $\bm{y}_{1}$ and the location of the sensors $\bm{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(\bm{y}_{2}|\bm{y}_{1},\bm{s})=p(\bm{y}_{2}|\bm{s})$.