Optimal Design of Bacterial Carpets for Fluid Pumping

Abstract

In this work, we outline a methodology for determining optimal helical flagella placement and phase shift that maximize fluid pumping through a rectangular flow meter above a simulated bacterial carpet. This method uses a Genetic Algorithm (GA) combined with a gradient-based method, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, to solve the optimization problem and the Method of Regularized Stokeslets (MRS) to simulate the fluid flow. This method is able to produce placements and phase shifts for small carpets and could be adapted for implementation in larger carpets and various fluid tasks. Our results show that given identical helices, optimal pumping configurations are influenced by the size of the flow meter. We also show that intuitive designs, such as uniform placement, do not always lead to a high-performance carpet.

1. Introduction

Bacterial carpets consist of multiple flagellated bacteria naturally or artificially adhered to a surface, with the flagella positioned outward into the fluid. The geometry of the surrounding fluid may be open [1] or it may be inside a tube or microscale channel [2]. The study of fluid flow or pumping by bacterial carpets has been of significant interest in microfluidics for the past couple of decades. There are numerous experimental studies on flow and mixing induced by helices [3,4,5], biomimetic cilia [6,7], and bacterial carpets [1,2,8,9,10].

These carpets can be modeled using various numerical representations of flagella by rigid or flexible helices. Because of the small length scales involved, these fluids are often modeled as low-Reynolds-number, Newtonian fluids [11]. Additionally, a number of numerical studies have been performed to examine the effect of these carpets on mixing and pumping [12,13,14,15]. Several of these studies have observed that helical placement and phase shifts alter the qualitative and quantitative characteristics of the fluid flow. However, to the best of the authors’ knowledge, no existing work has systematically looked at optimizing flow subject to these helical parameters.

Ideally these bacterial carpets are created in order to control or achieve either a directional flow or mixing in the fluid. In this study, the authors outline an approach to maximize the volume of fluid pumped through a flow meter placed slightly above the simulated helices subject to several variables such as the locations of the helices and their phase shifts. Such an optimization problem faces the following two fundamental challenges. First, a priori knowledge of what makes a carpet ‘good’ at pumping is lacking. Second, function evaluations for the objective task (pumping) are computationally costly fluid simulations. We seek a method that strikes a balance between reliability in finding the optimal carpet and computational efficiency.

The fluid simulations outlined in Section 2.1, Section 2.2 and Section 2.3 are based on those from [12,13] for simulating fluid flows induced by one or more rigid helical structures rotating with each helix attached to a wall (semi-infinite domain). Section 2.4 outlines a method which first employs a genetic algorithm [16,17,18] applied to a surrogate model in order to obtain a good starting guess for an optimal carpet. Compared to the full model, the surrogate model can be simulated at a substantially reduced computational cost. This carpet is then used as the initial solution in a gradient-based method [19,20,21,22,23,24,25] to find the optimal carpet. Our method takes advantage of the genetic algorithm’s ability to perform a global search along with the efficiency and accuracy of the gradient-based method. The procedure outlined here addresses optimizing the flow subject to helical parameters and has broad utility as a framework which may be adapted for a multitude of future applications. This method is tested and validated in Section 3 for cases of one, two, three, and four helices.

2. Methods

We model a bacterial carpet as a collection of rigid, helical filaments attached to a stationary planar wall and immersed in a viscous fluid. In Section 2.1, we lay out the Method of Regularized Stokeslets (MRS) used to solve the fluid equations. Section 2.2 and Section 2.3 describe the setup of both the rigid helices and the flow meter developed and utilized as a measure for the optimization described in Section 2.4.

2.1. Fluid Solver

The methods described in this section are pulled directly from [12,13]. The fluid is governed by the incompressible Stokes equations

$$\begin{array}{ccc}\hfill -\nabla p+\mu \Delta \mathbf{u}& =& \mathbf{F},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \nabla ·\mathbf{u}& =& 0,\hfill \end{array}$$

(2)

where $p\left(\mathbf{x}\right)$ is fluid pressure, $\mu$ is fluid viscosity, $\mathbf{u}\left(\mathbf{x}\right)$ is fluid velocity, and $\mathbf{F}\left(\mathbf{x}\right)$ is an external force density acting on the fluid. Using the MRS [26,27], we characterize the external force at each point ${\mathbf{x}}_k$ as a smooth concentrated force ${\mathbf{F}}_k\left(\mathbf{x}\right)={\mathbf{f}}_k{\psi }_{ϵ}\left(\right|\mathbf{x}-{\mathbf{x}}_k\left|\right)$. Here, ${\mathbf{f}}_k$ is a constant vector and ${\psi }_{ϵ}$ is a smooth approximation of the Dirac delta function commonly referred to as a blob function. In this work,

$${\psi }_{ϵ}\left(r\right)=\frac{15{ϵ}^4}{8\pi {(r^2+{ϵ}^2)}^{7/2}}.$$

The forces acting on the fluid come from the N points discretizing the rigid bodies in the fluid. Note that the MRS calculates the exact solution due to a regularized point force. To model a helical filament with nonzero effective thickness, we distribute regularized point forces along its centerline. Implementing the MRS, the velocity $\mathbf{u}$ at any point $\mathbf{x}$ resulting from forces ${\mathbf{f}}_k$ located at positions ${\mathbf{x}}_k$ is

$${\displaystyle \mathbf{u}\left(\mathbf{x}\right)=\sum _{k=1}^N\left[\frac{{\mathbf{f}}_k(r_k^2+2{ϵ}^2)+({\mathbf{f}}_k·(\mathbf{x}-{\mathbf{x}}_k))(\mathbf{x}-{\mathbf{x}}_k)}{8\pi {(r_k^2+{ϵ}^2)}^{3/2}}\right],}$$

(3)

where $r_k=|\mathbf{x}-{\mathbf{x}}_k|$. Note that the forces are implemented at the positions ${\mathbf{x}}_k$, and the velocity can be evaluated at any point $\mathbf{x}$ without requiring an underlying computational grid for the fluid domain, hence the velocity is defined everywhere.

When evaluating the velocity at the positions ${\mathbf{x}}_k$, Equation (3) can be written as

$$\mathcal{U}=A\mathcal{F}$$

(4)

where A is $3N\times 3N$, $\mathcal{U}$ is $3N\times 1$ and contains the N velocities $\mathbf{u}\left({\mathbf{x}}_k\right)$, and $\mathcal{F}$ is $3N\times 1$ and contains the N forces ${\mathbf{f}}_k$. To prescribe the kinematics of the rigid, rotating helices, one must define the velocities at the N points on the helices and solve the linear system (4) for the N forces that result in the desired velocity. These forces can then be used in Equation (3) to find the resulting velocity at any point in the unbounded domain. In practice, we solve a non-dimensional system, as is shown in detail in [12].

To incorporate a stationary planar wall, the method of images for regularized Stokeslets [28] is used to enforce a no-slip $(\mathbf{u}=\mathbf{0})$ boundary condition. In this formulation, an image source consisting of a regularized Stokeslet, doublet, dipole, and rotlet is mirrored across the boundary. This results in an infinite planar wall at which the velocity is zero, and no underlying grid is required to model the wall. In previous work [13], it was shown that the choice of boundary condition affects the efficiency of a bacterial carpet at pumping fluid. Therefore, we expect that when a different boundary condition is considered, such as the periodic boundary condition [29,30], the optimal configuration of the carpet will change. In terms of implementation, the optimization methods described in Section 2.4 still apply; the only change is how the objective function of the optimization problem (see Section 2.3) should be evaluated. (To be more precise, using a different boundary condition only affects how the flow field is calculated within the objective function.)

2.2. Helical Flagellum

As in previous work [12,13], we model each bacterial flagellum as a rotating helical filament, ${\mathbf{x}}_c=(x_c,y_c,z_c)$, whose centerline, parametrized by $\sigma$, is

$$\begin{array}{ccc}\hfill x_c\left(\sigma \right)& =& \alpha tanh\left(\tau \sigma \right)cos\left(\frac{2\pi \sigma }{\lambda }+\varphi \right),\hfill \\ \hfill y_c\left(\sigma \right)& =& \alpha tanh\left(\tau \sigma \right)sin\left(\frac{2\pi \sigma }{\lambda }+\varphi \right),\hfill \\ \hfill z_c\left(\sigma \right)& =& \sigma ,\hfill \end{array}$$

(5)

for $0\le \sigma \le L$. Here, L is the height of the helix, $\alpha$ is the helical radius, $\tau$ is a tapering parameter, $\lambda$ is the helical pitch, and $\varphi$ is the phase shift. We then discretize each helix into n equally spaced points. As in the Refs. [13,31], we choose parameters to be close to the typical dimensions of Escherichia coli bacteria. All default parameters, given in non-dimensional units, are listed in Table 1.

The kinematics of each helix are prescribed using the time-dependent rotation matrix

$$R_{\omega }=\left[\begin{array}{ccc}cos\left(\omega t\right)& -sin\left(\omega t\right)& 0\\ sin\left(\omega t\right)& cos\left(\omega t\right)& 0\\ 0& 0& 1\end{array}\right]$$

with angular velocity $(0,0,\omega )$. A negative value of $\omega$ results in a clockwise rotation when viewed from above. At any point in time, the position of the helix is given by

$$\begin{array}{c}\hfill \mathbf{x}\left(t\right)=R_{\omega }{\mathbf{x}}_c+{\mathbf{x}}_b,\end{array}$$

(6)

where ${\mathbf{x}}_b$ is the location of the base point.

In this work, the helices are assumed to be rigid. Under this assumption, we can prescribe their shape and motion as above. In many biophysical problems, the shape and motion of cilia or flagella are a result of fluid-structure interactions and typically unknown a priori. Models that are more suitable for these problems can be found in [32,33], for example. Although the model described here is not realistic in some applications, such as carpets made up of live bacteria, it is appropriate for man-made carpets and allows for a straight-forward investigation into questions of optimization.

2.3. Measurement of Pumping

In order to quantify fluid pumping by a carpet of helices described in Section 2.2, we measure the volume of fluid pumped upward through a flow meter stationed above the carpet during one rotation of the helices, as in the Ref. [13].

More precisely, the flow meter is chosen to be the following square “patch” parallel to the no-slip wall at $z=0$:

$$\left\{(x,y,z_m)|0\le x\le W\mathrm{and}0\le y\le W\right\},$$

(7)

where $z_m$ is greater than the height L of the helices and W denotes the width of the patch. Figure 1 shows such a flow meter above a carpet consisting of two helices. We note that the only boundary of the fluid is the infinite, planar, and solid wall at $z=0$, on which the flow vanishes. The flow meter simply provides a means of quantifying pumping and is not physically placed in the fluid. Thus, it has no effect on the fluid flow.

Let $\mathbf{\theta }\in {\mathbb{R}}^d$ contain all the design parameters that we aim to optimize. As specified in Section 3, in this work, the design parameters are the $(x,y)$ coordinates of the base points and the phase shifts of the helices. Nevertheless, the methods described in this section and in Section 2.4 still apply if different design parameters are considered. The following function then quantifies pumping for a carpet whose configuration is determined by the choice of $\mathbf{\theta }$:

$$f\left(\mathbf{\theta }\right)={\int }_0^W{\int }_0^W{\int }_0^T{u}_3\left(x,y,z_m,t,\mathbf{\theta }\right)dtdxdy,$$

(8)

where T is the time that it takes for the helices to perform one rotation and $u_3\left(x,y,z_m,t,\mathbf{\theta }\right)$ denotes the z- direction fluid velocity at a point $(x,y,z_m)$ on the flow meter (7) at time t. Given $\mathbf{\theta }$, we estimate $f\left(\mathbf{\theta }\right)$ as follows. We first approximate the triple integral in Equation (8) by the following double integral:

$${\int }_0^W{\int }_0^W\sum _{k=1}^{N_t}{u}_3\left(x,y,z_m,t_{k-1},\mathbf{\theta }\right)·\Delta tdxdy,$$

(9)

where $N_t$ is the number of time steps, $\Delta t=T/N_t$ is the size of the time step, and $t_{k-1}=(k-1)·\Delta t$ for $k=1,2,\dots ,N_t$; and then we evaluate (9) using a quadrature rule. The values of T and $N_t$ used in the numerical experiments can be found in Table 1. The fluid velocities $u_3$ at the quadrature points discretizing the flow meter (7) are calculated using the MRS.

As outlined in Section 2.1, at each time point $t_{k-1}$, calculating the forces exerted by the points on the helices to the surrounding fluid requires solving a linear system (4) whose coefficient matrix is of size $3N\times 3N$, where N is the number of points. The computational cost, using the standard “big $\mathcal{O}$” notation, is therefore $\mathcal{O}\left(N^3\right)$. After the forces have been calculated, evaluating $u_3\left(x,y,z_m,t_{k-1},\mathbf{\theta }\right)$ at $N_q$ quadrature points entails $\mathcal{O}\left(N_{q}N\right)$ additional operations (see Equation (3)). In summary, the computational complexity of evaluating $f\left(\mathbf{\theta }\right)$ is

$$\left(\mathcal{O}\left(N^3\right)+\mathcal{O}\left(N_{q}N\right)\right)·N_t.$$

(10)

This implies that the computational cost of evaluating $f\left(\mathbf{\theta }\right)$ increases cubically with the number of helices in the carpet, assuming that the same number of points are used to discretize each helix. The computational cost also increases with the size of the flow meter since more quadrature points are needed to approximate the integral (9) if the same accuracy is desired.

2.4. Optimization

To optimize the helical parameters of a bacterial carpet for a task (fluid pumping in the case of this work), we must solve the following optimization problem:

$$\mathrm{Find}{\mathbf{\theta }}^{*}\mathrm{that}\mathrm{maximizes}f\left(\mathbf{\theta }\right),$$

(11)

where $\mathbf{\theta }\in {\mathbb{R}}^d$ contains d design parameters, such as the locations of the helices and their phase shifts, $f\left(\mathbf{\theta }\right):{\mathbb{R}}^d\to \mathbb{R}$ is a scalar function that measures the performance of the carpet at the task. (Depending on the task and definition of $f\left(\mathbf{\theta }\right)$, we may wish to maximize or minimize $f\left(\mathbf{\theta }\right)$. In this work, the goal is to maximize $f\left(\mathbf{\theta }\right)$). Unfortunately, an analytic form of $f\left(\mathbf{\theta }\right)$ is usually unavailable and evaluating it entails numerical simulations of the fluid around the carpet. To solve (11), we need to apply a numerical optimization algorithm, which in turn requires evaluating f at different $\mathbf{\theta }$’s. In this section, we discuss two such algorithms: Genetic Algorithms (GAs) which are gradient-free and typically used to find a global optimizer, and the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, a gradient-dependent, quasi-Newton method that can converge to a local optimizer quickly if a good estimate is given. To find a global optimizer of $f\left(\mathbf{\theta }\right)$ efficiently, we propose to use a hybrid method that exploits both a GA’s ability at identifying a global optimizer and the BFGS algorithm’s fast convergence.

2.4.1. Genetic Algorithms (GAs)

GAs [16,17,18] are often applied to find the global minimum or maximum of a function. They “evolve” a pool of candidate solutions towards better solutions through an iterative process inspired by natural selection. More specifically, the best solutions of the current generation are selected, recombined, and possibly mutated to produce the next generation of solutions. We first randomly select the initial generation of solutions. In each iteration that follows, a portion of the current generation of solutions is chosen at random to be “parents” to the next generation. The solutions that produce lower (higher) function values are more likely to be chosen, if the minimum (maximum) is sought. Two parent solutions can be combined as follows to produce a “child” solution. Every entry of the child solution is set to the corresponding entry of one of the parents, depending on the result of a coin flip. We can also mutate a parent solution to produce a child solution by adding a number taken from a Gaussian distribution with mean 0 to each entry of the parent solution. GAs return the best solution among the final generation of candidate solutions as the optimizer. We note that the methods described above are not unique, and there exist many options for the creation of the initial population and the selection, crossover, and mutation of the parent solutions. Custom methods that are tailored to the optimization problem (11) at hand can also be devised. Like natural selection, GAs are stochastic, and therefore, their performance may be affected by the sequence of random numbers used.

GAs and gradient-free optimization methods in general suffer from slow convergence. As will be shown in Section 3.1, even for the simplest carpet made up of a single helix, a GA entails thousands of evaluations of the objective function. This can be prohibitively expensive for fluid problems such as the ones considered here.

2.4.2. The Broyden-Fletcher-Goldfarb-Shanno (BFGS) Algorithm

The BFGS algorithm [21,22,23,24,25] is an iterative method for optimization that starts with an initial guess ${\mathbf{\theta }}_0$ and calculates a new estimate as

$${\mathbf{\theta }}_{k+1}={\mathbf{\theta }}_k+{\alpha }_k{\mathbf{\delta }}_k,for k=0,1,2,\dots ,$$

(12)

where ${\alpha }_k$ is a scalar and ${\mathbf{\delta }}_k$ a search direction. Like in Newton’s method, the search direction ${\mathbf{\delta }}_k$ in Equation (12) is obtained by solving a linear system

$$H_k{\mathbf{\delta }}_k=\nabla f\left({\mathbf{\theta }}_k\right),$$

(13)

where $\nabla f\left({\mathbf{\theta }}_k\right)$ is the gradient of f at ${\mathbf{\theta }}_k$. In Newton’s method, the coefficient matrix $H_k$ in Equation (13) is the Hessian matrix associated with f at ${\mathbf{\theta }}_k$, whereas the BFGS algorithm uses an approximation to the Hessian instead that is considerably easier to calculate and invert. After ${\mathbf{\delta }}_k$ has been calculated, ${\alpha }_k$ in Equation (12) is then chosen to be the scalar that optimizes $f\left({\mathbf{\theta }}_k+\alpha {\mathbf{\delta }}_k\right)$. When an analytic form of $f\left(\mathbf{\theta }\right)$ is unavailable, such as in the problems considered here, the gradient $\nabla f\left({\mathbf{\theta }}_k\right)$ needs to be approximated numerically, for example, by the forward difference method.

The BFGS algorithm converges rapidly to a local optimizer of $f\left(\mathbf{\theta }\right)$ if the initial guess ${\mathbf{\theta }}_0$ is close to it. However, the algorithm may get “stuck” at a local optimizer that is not global if ${\mathbf{\theta }}_0$ is not carefully chosen; unfortunately, for the problems considered here, it remains unclear how to choose it. As will be shown in Section 3, even for simple carpets, the initial configuration that we intuit may turn out to be inferior.

2.4.3. The Hybrid Method

To take advantage of both a GA’s higher reliability at finding the global optimizer and the BFGS algorithm’s higher efficiency given a good initial guess, we propose to use a hybrid approach for solving (11) that consists of the following two steps:

Step 1:

find a maximizer ${\mathbf{\theta }}^{†}$ of $s\left(\mathbf{\theta }\right)$ using a GA, where s is a surrogate model for f;

Step 2:

find a maximizer ${\mathbf{\theta }}^{*}$ of $f\left(\mathbf{\theta }\right)$ using the BFGS algorithm with ${\mathbf{\theta }}^{†}$ as the initial guess for ${\mathbf{\theta }}^{*}$.

The surrogate model $s\left(\mathbf{\theta }\right)$ should be significantly cheaper to evaluate than $f\left(\mathbf{\theta }\right)$ so that it is feasible to apply a GA to it in Step 1, yet it should also be sufficiently representative of $f\left(\mathbf{\theta }\right)$ so that ${\mathbf{\theta }}^{†}$ will be close to ${\mathbf{\theta }}^{*}$ and thus the BFGS algorithm will converge quickly in Step 2. The $s\left(\mathbf{\theta }\right)$ that we use employs coarser temporal, spatial discretization and will be specified in Section 3.

3. Results

All the computational results are obtained in MATLAB^® version R2017a on a virtual machine equipped with two Intel^® Xeon^® E7-8890 v4 Processors with twenty-four cores per processor and 64GB of RAM. MATLAB’s Parallel Computing Toolbox is used to parallelize the computations. Default parameters are given in

Table 1. Default parameter values given in non-dimensional units.

Parameter	Description	Value
$ϵ$	blob size	0.01 (for f) or 0.1 (for s)
$\mu$	viscosity	1
L	axial length of helix	2.2
$\alpha$	helical radius	0.085
$\tau$	tapering parameter	1000
$\lambda$	helical pitch	$2\pi \alpha$
$\omega$	angular speed	$-2\pi$
n	points per helix	158 (for f) or 16 (for s)
T	time to perform one rotation	$2\pi /\left|\omega \right|$
$N_t$	number of time steps per rotation	1000 (for f) or 100 (for s)
$z_m$	flow meter location	2.5

and do not vary throughout the numerical experiments. Using a different set of values will only affect the evaluation of the objective function and does not affect how the optimization methods described in Section 2.4 will be implemented. The design parameters are chosen to be the phase shifts and $(x,y)$ coordinates of the bases of the helices, and we aim to determine their optimal values that achieve maximal pumping using the methods in Section 2.4. We consider carpets made up of a single helix (Section 3.1), a pair of helices (Section 3.2), and multiple helices (Section 3.3).

Applying the hybrid method described in Section 2.4.3 to the optimization problem (11) entails evaluating the objective function $f\left(\mathbf{\theta }\right)$ and its surrogate $s\left(\mathbf{\theta }\right)$ repeatedly at different $\mathbf{\theta }$’s. As shown in Table 1, the spatial and temporal discretization used in s is about ten times as coarse as the discretization used in f. For flow meters of various sizes and carpets consisting of various numbers of helices, we compare the runtimes required to evaluate $f\left(\mathbf{\theta }\right)$ and $s\left(\mathbf{\theta }\right)$. Specifically, we consider three values of the width W of the flow meter (7): $2.5$, 5, and 10 as well as cases of one, two, and four helices. The runtimes entailed by evaluating $f\left(\mathbf{\theta }\right)$ and $s\left(\mathbf{\theta }\right)$ for each combination of W and number of helices are listed in

Table 2. The runtimes of evaluating $f\left(\mathbf{\theta }\right)$ and $s\left(\mathbf{\theta }\right)$ when the carpet consists of one, two, or four helices and the width of the flow meter (7) is $2.5$, 5, or 10.

W	One Helix		Two Helices		Four Helices
	$\mathit{f}$	$\mathit{s}$	$\mathit{f}$	$\mathit{s}$	$\mathit{f}$	$\mathit{s}$
$2.5$	165	$9.8$	449	$17.1$	763	$25.3$
5	244	$13.4$	823	$23.9$	1365	$44.8$
10	489	$19.4$	1416	$40.3$	2494	$79.8$

. Here, and for the rest of the paper, runtimes are measured in seconds by MATLAB’s tic and toc commands. Recall that a double integral must be evaluated to calculate $f\left(\mathbf{\theta }\right)$ or $s\left(\mathbf{\theta }\right)$ (see (9)). We use MATLAB’s function integrate2 with its default setting to approximate this integral numerically.

As seen in Table 2, the time to evaluate both $f\left(\mathbf{\theta }\right)$ and $s\left(\mathbf{\theta }\right)$ increases with the number of helices and the size of the flow meter; this is expected given the computational complexity (10). The runtime increases with the size of the flow meter because more quadrature points are needed as the domain of integration in (9) increases. For the same flow meter and number of helices, evaluating $s\left(\mathbf{\theta }\right)$ is significantly cheaper than evaluating $f\left(\mathbf{\theta }\right)$.

3.1. A Single Helix

We first consider the case of a single helix by investigating the following question: Where should the helix be placed to maximize pumping? We fix all design parameters of the helix except for the $(x,y)$ coordinates of its base point, that is, $\mathbf{\theta }=\left(x_{\mathrm{b}},y_{\mathrm{b}}\right)\in {\mathbb{R}}^2$ in the optimization problem (11).

In Figure 2, we show the streamlines drawn based on the instantaneous fluid velocity around a rotating helix. As can be seen in Figure 2d,e, the helix moves the fluid near its tip upward while rotating, and there are also regions of downward flow away from the helix. A more detailed experimental analysis of the flow field around a single rotating helix can be found in [31] and a numerical study was performed in [12,13].

Since we consider an incompressible Newtonian fluid, as the width W of the flow meter (7) approaches infinity, the volume of the fluid pumped through the meter goes to zero regardless of the values of the design parameters. However, in applications, only flow meters of finite sizes are of interest, and the volume of fluid that is transported through a finite flow meter does change depending on the placement of the helix relative to the flow meter.

We consider a square flow meter defined by (7) and three values of W: $2.5$, 5, and 10. In Figure 3b,d,f, we show the surface and contour plots of $f\left(\mathbf{\theta }\right)$ over the domain $[0,W]\times [0,W]$ corresponding to $W=2.5$, 5, and 10, respectively. For each W, in order to produce the pair of plots, we sample f at the following uniformly distributed grid points:

$${\left\{\left((i-1)·\frac{W}{100},(j-1)·\frac{W}{100}\right)\right\}}_{i,j=1}^{101}.$$

(14)

As can be seen in Figure 3b,d,f, for all three values of W, global maximizers of f lie on the diagonals of the domain. It is interesting to note that while the global maximizer is the center $\left(\frac{W}{2},\frac{W}{2}\right)$ of the domain and unique when $W=2.5$ or 5, there are four global maximizers located near the four corners of the domain when $W=10$. While intuition might lead one to expect the center to be the global maximizer, in fact, it is a local minimizer when $W=10$, and the global maximum ($0.0389$) is considerably larger than this local minimum ($0.0250$). This clearly demonstrates the importance of quantitative analysis when searching for an optimal design of a bacterial carpet.

To further investigate why the optimal location of the helix is near the corners instead of the center of the domain when the width of the flow meter is 10, in Figure 4, we draw the contour plots of the instantaneous z-direction fluid velocity on the flow meter when the helix is based at one of the four optimal locations (Figure 4a) and when it is based at the center of the domain (Figure 4b). Comparing the two contour plots, we observe that when the helix is based at the center, a considerably larger region on the flow meter has small negative z-direction velocity. More precisely, the areas enclosed by the level curves corresponding to z-direction velocities $-0.00025$, $-0.0005$, $-0.00075$, and $-0.0009$ in Figure 4b are almost twice as large as their counterparts in Figure 4a. In addition, Figure 4 also helps explain why the optimal location of the helix is at the center when the width of the flow meter is $2.5$ or 5: both flow meters are small enough that placing the helix at the center maximizes the region with large positive z-direction velocity while creating no region or only a negligible region where the z-direction velocity is negative.

We apply both the hybrid method described in Section 2.4.3 and a GA to solve (11) and compare their performance. In Step 1 of the hybrid method, we apply a GA to find a maximizer ${\mathbf{\theta }}^{†}$ of $s\left(\mathbf{\theta }\right)$, the surrogate model of $f\left(\mathbf{\theta }\right)$ (see Table 1 and Table 2 for a comparison of the two functions); and in Step 2 of the hybrid method, we apply the BFGS algorithm with initial guess ${\mathbf{\theta }}^{†}$ to find a maximizer ${\mathbf{\theta }}^{*}$ of $f\left(\mathbf{\theta }\right)$. In

Table 3. The maximizers found by the hybrid method for the case of a single helix and three sizes of the flow meter (7): $W=2.5$, 5, and 10.

W	Step 1: GA Applied to s	Step 2: BFGS Applied to f
	$\left({\mathbf{x}}_{\mathbf{b}}^{†},{\mathbf{y}}_{\mathbf{b}}^{†}\right)$	$\left({\mathbf{x}}_{\mathbf{b}}^{*},{\mathbf{y}}_{\mathbf{b}}^{*}\right)$
$2.5$	$(1.2679,1.2630)$	$(1.2500,1.2500)$
5	$(2.3757,2.5106)$	$(2.5000,2.5000)$
10	$(1.7245,8.2893)$	$(1.7104,8.2895)$

, we list the maximizers ${\mathbf{\theta }}^{†}$ and ${\mathbf{\theta }}^{*}$ found in the hybrid method for flow meters (7) of widths $2.5$, 5, and 10. In Figure 3a,c,e, we also display the top views of the optimal carpets corresponding to the three flow meters respectively.

Table 4. The runtimes, iteration counts, and numbers of function evaluations entailed by the hybrid method for the case of a single helix and three sizes of the flow meter (7): $W=2.5$, 5, and 10. (Forty-eight computer cores are used).

W	Step 1: GA Applied to s			Step 2: BFGS Applied to f			Total
	# of Iter.	# of s Eval.	Runtime	# of Iter.	# of f Eval.	Runtime	Runtime
$2.5$	66	3350	877	2	12	207	1084
5	53	2700	1170	3	18	386	1556
10	64	3250	1717	2	15	461	2178

and

Table 5. The maximizers found by the GA and the runtimes, iteration counts, and numbers of function evaluations entailed by the GA for the case of a single helix and three sizes of the flow meter (7): $W=2.5$, 5, and 10. (Forty-eight computer cores are used).

W	$\left({\mathit{x}}_{\mathit{b}}^{**},{\mathit{y}}_{\mathit{b}}^{**}\right)$	# of Iter.	# of f Eval.	Runtime
$2.5$	$(1.2679,1.2630)$	67	3400	19,639
5	$(2.4461,2.4573)$	55	2800	24,784
10	$(1.7245,1.6907)$	64	3250	38,741

contain the numbers of iterations, the numbers of function evaluations, and the runtimes of the hybrid method and a GA, respectively. We use MATLAB’s ga function, an implementation of a GA, to find maximizers for both $s\left(\mathbf{\theta }\right)$ and $f\left(\mathbf{\theta }\right)$. The default setting of ga is used except that parallel computing is enabled and the range of the initial generation of candidate solutions is set to $[0,W]\times [0,W]$. We note that GAs are not deterministic (see Section 2.4.1). Thus, for reproducibility of the results presented in Table 3, Table 4 and Table 5, we reset the seed of MATLAB’s random number generators to 0 using the command rng default before calling the function ga. We use MATLAB’s fminunc function with its default setting, which implements the BFGS algorithm, to find a maximizer of $f\left(\mathbf{\theta }\right)$ in Step 2 of the hybrid method (Both MATLAB functions ga and fminunc are for the minimization of a given objective function. Since we wish to maximize $s\left(\mathbf{\theta }\right)$ or $f\left(\mathbf{\theta }\right)$ instead, in our MATLAB implementation, we use $-s\left(\mathbf{\theta }\right)$ or $-f\left(\mathbf{\theta }\right)$ as the objective function).

As seen in Table 4 and Table 5, for all three flow meters, the hybrid method is significantly more efficient than the GA, requiring $90\%$ less runtime. Although the GA requires about three thousand function evaluations to converge whether it is applied to $s\left(\mathbf{\theta }\right)$ or $f\left(\mathbf{\theta }\right)$, since $s\left(\mathbf{\theta }\right)$ is considerably cheaper to evaluate (see Table 2), applying the GA to $s\left(\mathbf{\theta }\right)$ is also considerably cheaper. In the hybrid method, once the estimate ${\mathbf{\theta }}^{†}$ is found by the GA in Step 1, the BFGS algorithm, which starts with ${\mathbf{\theta }}^{†}$, converges quickly to the maximizer ${\mathbf{\theta }}^{*}$ in Step 2, requiring only two or three iterations and less than 20 function evaluations. Moreover, we can see that the hybrid method is more accurate than the GA by comparing the maximizers $\left(x_b^{*},y_b^{*}\right)$ and $\left(x_b^{**},y_b^{**}\right)$ found by the two methods, respectively (see Table 3 and Table 5). We also observe that the runtimes of both the hybrid method and GA increase with the size of the flow meter. Note that since the initial candidate solutions are randomly chosen in the GA, in the case where the flow meter is $10\times 10$, both the hybrid method and GA may find different maximizers in different runs if we do not reset the seed of the random number generators to the same number.

To demonstrate the robustness of the GA at finding a good initial guess in Step 1 of the hybrid method, for the case where the flow meter is $10\times 10$, we apply the GA to $s\left(\mathbf{\theta }\right)$ sixty-four times and initialize the random number generators based on the current time using the command rng shuffle before each run. The GA uses a different sequence of random numbers in each run due to the re-initialization. In Figure 5a, we display all sixty-four ${\mathbf{\theta }}^{†}$’s found by the GA. They form four clusters around the four maximizers of $f\left(\mathbf{\theta }\right)$ (see Figure 3f), and the sizes of the clusters are 20, 19, 10, and 15. This demonstrates that the GA applied to $s\left(\mathbf{\theta }\right)$ always finds an initial guess ${\mathbf{\theta }}^{†}$ close to one of the maximizers of $f\left(\mathbf{\theta }\right)$, although the maximizer may differ from run to run depending on the sequence of random numbers used. Runtime of the GA varies with the seed of the random number generators as well. In Figure 5b, we show a histogram of the sixty-four runtimes, the majority of which is less than 1800 s.

3.2. A Pair of Helices

We now apply the hybrid method to find an optimal design of a carpet consisting of two helices for fluid pumping. We set the phase shift of the first helix to zero (${\varphi }_1=0$) and let the design parameters be the $(x,y)$ coordinates of the base points of the two helices and the phase shift of the second helix. That is,

$$\mathbf{\theta }=\left(x_b^1,y_b^1,x_b^2,y_b^2,{\varphi }_2\right)\in {\mathbb{R}}^5.$$

(15)

This optimization problem (11) is significantly more difficult to solve than the one in Section 3.1, not only because there are more design parameters (five instead of two) but also because the computational cost of evaluating $f\left(\mathbf{\theta }\right)$ and $s\left(\mathbf{\theta }\right)$ increases with the number of helices (see (10)) as is shown in Table 2.

In Figure 6, we show the streamlines drawn based on the instantaneous fluid velocity around a pair of rotating helices. As can be seen in Figure 6d,e, the helices move the fluid near their tips upward while rotating. A more detailed analysis of the flow field around a pair of rotating helices can be found in [13].

We again consider a square flow meter (7) and three values of W, as in Section 3.1. The maximizers found by the hybrid method are listed in

Table 6. The maximizers found by the hybrid method for the case of a pair of helices and three sizes of the flow meter (7): $W=2.5$, 5, and 10.

W	Step 1: GA Applied to s			Step 2: BFGS Applied to f
	$\left({\mathbf{x}}_{\mathbf{b}}^{\mathbf{1}†},{\mathbf{y}}_{\mathbf{b}}^{\mathbf{1}†}\right)$	$\left({\mathbf{x}}_{\mathbf{b}}^{\mathbf{2}†},{\mathbf{y}}_{\mathbf{b}}^{\mathbf{2}†}\right)$	${\mathbf{\varphi }}_{\mathbf{2}}^{†}$	$\left({\mathbf{x}}_{\mathbf{b}}^{\mathbf{1}*},{\mathbf{y}}_{\mathbf{b}}^{\mathbf{1}*}\right)$	$\left({\mathbf{x}}_{\mathbf{b}}^{\mathbf{2}*},{\mathbf{y}}_{\mathbf{b}}^{\mathbf{2}*}\right)$	${\mathbf{\varphi }}_{\mathbf{2}}^{*}$
$2.5$	$(1.6099,0.8510)$	$(0.9640,1.6778)$	$3.0700$	$(1.5405,0.9276)$	$(0.9595,1.5724)$	$3.1420$
5	$(3.1618,1.7583)$	$(1.9989,2.9001)$	$3.2090$	$(3.0124,1.9583)$	$(1.9877,3.0416)$	$3.1410$
10	$(1.6746,8.3083)$	$(8.3767,1.8351)$	$3.1491$	$(1.7091,8.2913)$	$(8.2909,1.7087)$	$3.1491$

and the top views of the corresponding optimal carpets are displayed in Figure 7. We use ${\mathbf{\theta }}^{†}=\left(x_b^{1†},y_b^{1†},x_b^{2†},y_b^{2†},{\varphi }_2^{†}\right)$ to denote the maximizer of $s\left(\mathbf{\theta }\right)$ found by the GA and ${\mathbf{\theta }}^{*}=\left(x_b^{1*},y_b^{1*},x_b^{2*},y_b^{2*},{\varphi }_2^{*}\right)$ to denote the maximizer of $f\left(\mathbf{\theta }\right)$ found by the BFGS algorithm using ${\mathbf{\theta }}^{†}$ as its initial guess. The iteration counts, numbers of function evaluations, and runtimes of the GA and BFGS algorithm are reported in

Table 7. The runtimes, iteration counts, and numbers of function evaluations entailed by the hybrid method for the case of a pair of helices and three sizes of the flow meter (7): $W=2.5$, 5, and 10. (Forty-eight computer cores are used).

W	Step 1: GA Applied to s			Step 2: BFGS Applied to f			Total
	# of Iter.	# of s Eval.	Runtime	# of Iter.	# of f Eval.	Runtime	Runtime
$2.5$	105	5300	2174	19	132	2595	4769
5	68	3450	2021	17	114	4273	6294
10	120	6050	6229	4	42	2565	8794

.

As can be seen in Table 6 and Figure 7, the hybrid method suggests that in order to maximize pumping, we should place the pair of helices along the diagonal of the domain $[0,W]\times [0,W]$ and set the phase shift of the second helix to $\pi$ (recall that the phase shift of the first helix is fixed at 0). These results are consistent with what was observed in [13] in that pumping was maximized when the difference between the phase shifts is $\pi$. However, we note that while this does maximize pumping, the distance between the optimal locations of the helices is sufficiently large so that when they are placed at these locations, the choice of phase shifts has a minimal impact on pumping. (For a pair of helices, the effect of the difference between their phase shifts on pumping has been studied in [13], where they were placed at various distances apart.) This is similarly observed for the cases of three and four helices considered in Section 3.3. We conjecture that the choice of phase shifts will play a greater role when the optimal locations of the helices are forced to be closer to one another, such as when there are more helices or when the flow meter is smaller.

Table 7 shows that, like in the case of a single helix, although a large number of iterations and function evaluations are needed for the GA to find the initial guess, ${\mathbf{\theta }}^{†}$, once it is found, the BFGS algorithm only requires a relatively small number of iterations and function evaluations to converge to the maximizer, ${\mathbf{\theta }}^{*}$. Additionally, the total runtime of the hybrid method increases with the size of the flow meter.

In the case with one helix, because we were optimizing over only two parameters, the $(x,y)$ coordinates of the base of the helix, we could easily visualize the objective function. Now however, we are optimizing over five parameters (the base coordinates of each helix, and the phase shift of the second) so it is more difficult to visualize our objective function in a useful way. In order to investigate how the placement of the pair of helices affects the carpet’s performance at pumping, for the flow meters (7) with widths $W=5$ and 10, we fix the phase shifts of both helices, the location of one of the them and calculate $f\left(\mathbf{\theta }\right)$ as the location of the other helix is varied. To be more precise, we fix all the design parameters of the optimal carpets shown in Figure 7b,c except for the location of the lower helix. We then examine the value of $f\left(\mathbf{\theta }\right)$ as we move the lower helix along the dashed line shown in Figure 8a,c. This line goes through the base points $\left(x_b^{1*},y_b^{1*}\right)$ and $\left(x_b^{2*},y_b^{2*}\right)$ of the two helices in the optimal carpet found by the hybrid method. As can be seen in Figure 8b, when $W=5$, the pumping capability of the carpet first grows as the distance increases, peaks as the moving helix reaches its optimal position $\left(x_b^{1*},y_b^{1*}\right)$ found by the hybrid method, and then starts to decline as the distance increases further. For the case $W=2.5$, the graph of $f\left(\mathbf{\theta }\right)$ is quite similar and not shown here. As shown in Figure 8d, the graph of $f\left(\mathbf{\theta }\right)$ is even more interesting in the case of $W=10$. As the distance increases, the value of $f\left(\mathbf{\theta }\right)$ first increases to reach a local maximum and then decreases to reach a local minimum, before it increases again to the global maximum when the moving helix is at its optimal position $\left(x_b^{2*},y_b^{2*}\right)$ calculated by the hybrid method. In Figure 8, a few sample locations of the moving helix as well as the corresponding values of $f\left(\mathbf{\theta }\right)$ are marked for reference.

3.3. Multiple Helices

We also apply the hybrid method to determine optimal designs that maximize pumping for carpets consisting of three and four helices. The design parameters are the $(x,y)$ coordinates of the base points and phase shifts of the helices. (The phase shift of the first helix is fixed at 0, as in Section 3.2). Thus, there are eight and eleven design parameters in these two cases respectively.

In the top two rows and bottom two rows of Figure 9, we show the streamlines drawn based on the instantaneous fluid velocity around three and four rotating helices, respectively. As can be seen in Figure 9d,e,i,j, the helices move the fluid near their tips upward while rotating.

We again consider the square flow meter (7) whose width is $W=5$ and display the top views of the optimal carpets in Figure 10. Interestingly, the helices approximately form an equilateral triangle when there are three of them (Figure 10a) and a square when there are four (Figure 10b) and both shapes are centered at $\left(2.5,2.5\right)$. The differences between the phase shifts of two helices located on adjacent vertices are $\frac{2\pi }{3}$ and $\pi$ in the two cases respectively.

For the case of four helices, we also examine the objective function $f\left(\mathbf{\theta }\right)$ for a family of carpets which includes the optimal one found by the hybrid method. They share the following properties: the base points of the four helices form a square centered at $\left(2.5,2.5\right)$, each of its four edges is parallel to the x or y axis, and the difference in the phase shifts of two helices located on adjacent vertices is $\pi$. In Figure 11a, we illustrate three such carpets whose base points are marked by crosses. The optimal carpet found by the hybrid method is also plotted for reference. For this family of carpets, we plot the value of $f\left(\mathbf{\theta }\right)$ against the width of the square in Figure 11b, which indicates that as the square becomes larger, the ability of the carpet to pump first increases and then decreases, reaching its maximum at the width predicted by the hybrid method. One key observation can be made from Figure 11 that when the four helices are uniformly distributed within the domain $[0,5]\times [0,5]$, that is, when the width of the square is $2.5$, the pumping is sub-optimal. Once again, this shows that intuition does not necessarily lead to optimal carpet design.

In

Table 8. The runtimes, iteration counts, and numbers of function evaluations entailed by the hybrid method for the cases of three and four helices. The flow meter is (7) with $W=5$. (Forty-eight computer cores are used).

# of	Step 1: GA Applied to s			Step 2: BFGS Applied to f			Total
Helices	# of Iter.	# of s Eval.	Runtime	# of Iter.	# of f Eval.	Runtime	Runtime
three	78	15,800	12,479	6	136	8615	21,094
four	159	32,000	30,178	38	450	45,729	75,907

, we report the numbers of iterations, numbers of function evaluations as well as runtimes of the hybrid method. Comparing Table 4, Table 7 and Table 8, we note that the computational cost of the hybrid method increases progressively with the number of helices. This is again due to the growing number of design parameters and cost of evaluating $f\left(\mathbf{\theta }\right)$ and $s\left(\mathbf{\theta }\right)$ (see (10) and Table 2).

4. Discussion

We outlined and validated a numerical method for optimizing a fabricated bacterial carpet for fluid pumping through a rectangular flow meter placed just above the helix or helices representing the bacterial flagella. The optimal placement of bacteria in a carpet has many practical applications in biofluidics, as well as the study of food capture or transport of nutrients for or by microorganisms [1,2,8,9,34,35]. To the authors’ knowledge, this is the first study to investigate optimization of carpet design for a flow task or objective. While other methods have previously numerically studied the effects of relative helical positions and phase shifts on pumping [13], this method is the first to study the selection of helical placement and phase shift to maximize fluid pumping.

The optimization problems considered here are challenging for several reasons. Firstly, local maximizers that are not global do exist, and gradient-based methods can converge quickly to these, missing a global maximizer. Secondly, the chosen objective function is not explicitly known, and evaluating it entails computationally expensive fluid simulations. Our approach first employs a GA to perform a global search for an initial solution, then uses this as the starting point for a gradient-based method to perform local searches for the optimal carpet. The GA is applied to a surrogate model for the objective function, which uses coarser temporal and spatial discretization to reduce computational costs. Once an initial solution is found by the GA, the gradient-based method is able to quickly and accurately identify a maximizer in its local vicinity.

The framework developed here is versatile and not limited to optimizing pumping or by the measurement tool employed here (the flow meter). It could also be adapted to consider more design parameters. That is, the carpet could be optimized for mixing or any other fluid task for which an objective function can be determined, or for other helical parameters. For example, since tilt plays an important role in pumping with synthetic flagella for biomedical devices [12,36], we might include tilt in future studies. This would introduce new complications since, for example with multiple helices, solving the optimization problem will be more complicated to avoid flagella physically overlapping.

Our results show that, with the current measurement of pumping, the optimal carpet design depends significantly on the size of the flow meter. Intuitive designs, such as uniform placement, do not always lead to a high-performance carpet. These results highlight the importance of quantitative analysis.

We observed that computing an optimal carpet becomes substantially more expensive as the number of helices increases. At one level, the efficiency of the GA could be improved by using custom methods for crossover and mutation that are more tailored to bacterial carpets. At another level, we want to reduce the computational costs associated with fluid simulation. For example, we might consider implementing regularized Stokeslet segments to reduce the number of points used to model the helices [37]. Additionally, when the number of helices is large, more sophisticated methods must be applied to solve for the forces exerted by the helices on the fluid [38] and calculate the fluid velocities from helical forces [39].