# Bayesian Learning Designs and Characterizes Porous Metamaterials for Biofilm
Transport and Control

Hanfeng Zhai, Jingjie Yeo
Sibley School of Mechanical and Aerospace Engineering
Cornell University, Ithaca, NY 14850, USA Corresponding Author:
<EMAIL_ADDRESS>
###### Abstract

Biofilm growth and transport in confined systems is a common phenomenon. While
machine learning (ML) and optimization have been extensively applied in
materials design, there is still a scarcity of thorough evaluations regarding
the optimization process. We combined Bayesian optimization (BO) and
individual-based modeling to conduct design optimizations for maximizing
different porous materials’ (PM) biofilm transporation capability. We first
characterize the acquisition function in BO for designing 2-dimensional porous
membranes. Results showed that the variance of the overall samples by the
upper confidence bound (UCB) is 32.08% higher than that of the expected
improvement (EI); the mean objective of the overall samples by the EI is 1.49%
higher than that of the UCB. Given the predefined target region, the EI is
2.35% more efficient than the UCB compared with uniform grid search. We then
use EI for designing lattice metamaterials (LM) and 3-dimensional porous media
(3DPM). It is found that BO is 92.89% more efficient than the uniform grid
search for LM and 223.04% more efficient for 3DPM. The selected
characterization simulation tests match well with the Gaussian process
regression approximated design spaces for three cases. We found that all the
extracted optimal designs exhibit better biofilm growth and transportability
than nonconfined vacuum space. Our comparison study shows that PM stimulates
biofilm growth by taking up volumetric space and pushing biofilms’ upward
growth, as evidenced by a 20% increase in biofilms in vacuum space compared to
porous materials. There are 128% more biofilms in the target growth region for
the PM-induced biofilm growth compared with the vacuum space growth. Our work
provides deeper insights into bio-porous materials design, ML and optimization
process characterizations, and extracting new physical mechanisms from the
optimizations.

_Keywords_ Individual-based modeling $\cdot$ biofilm $\cdot$ metamaterials
$\cdot$ Bayesian optimization $\cdot$ porous media

##  1 Introduction

Biofilms, commonly defined as surface-attached communities of microorganisms
(i.e., groups of bacteria cells) embedded in a self-produced matrix of
extracellular polymeric substances (EPS) Costerton et al. (1999), grow mostly
in confined systems such as rock cracks, industrial pipelines, biological
bodies, and many other artificial or natural microenvironments Friedlander et
al. (2013). One of the prerequisites of biofilm growth is the existence of
adhesive surfaces that allow bacteria to grow and clustered into “film-shaped”
communities adhered to by EPS. This accounts for the phenomenon that biofilms
are observed in mostly confined systems: an increasing surface area would
allow biofilm to attach and further grow. Biofilms, from the engineering
perspective, possess abundant pros and cons to human society. On the negative
side, the formation and attachment of biofilms pose serious problems for
marine engineering as cause fouling on the surfaces of marine vessels,
equipment, and infrastructure, leading to reduced efficiency and increased
maintenance costs Yebra et al. (2004); Dobretsov et al. (2006) and in
biomedical treatments as forming on medical devices, such as catheters and
implants, leading to infections that are difficult to treat Costerton et al.
(1999); Donlan and Costerton (2002). On the positive side, biofilm can also be
utilized as engineered living materials (ELM) to be used to create self-
healing concrete by incorporating bacteria into the concrete mix Jonkers et
al. (2010), treat wastewater by removing pollutants and nutrients
Chattopadhyay et al. (2022), and conduct 3D bioprinting for functional soft
materials Balasubramanian et al. (2019).

Considering all these pros and cons, one deduces that understanding the
mechanism of biofilm growth within confined systems for biofilm utilization is
significant for human-desired biofilm control. To summarize and elaborate on
this significance, there are three major points: (1) One can prevent undesired
biofilm attachment and conduct efficient biofilm removal Zhai and Yeo (2022).
(2) One may utilize the physics to promote efficient usage of biofilm as ELM,
e.g., clean energy applications Liu et al. (2022). (3) By combining both the
pros and cons to enable biofilm control to design customized devices and
sensors Mukherjee and Cao (2020). From both three points, we identify a major
goal that may potentially bring in the solutions for biofilm control: design
porous structural materials that can control biofilm growth. Following our
points, two problems naturally arise as efficient biofilm control and
utilization, (1) Conducting experiments on biofilm is time-consuming, making
it extremely difficult the characterize and benchmark in a timely effort; (2)
Directly modify the structures of porous materials to test the corresponding
biofilm growth properties is not straightforward and making the investigation
even more time-consuming. Hence, novel techniques that can bypass this “trial-
and-error” approach is of urgent need.

To tackle the first problem, our solution is to use computational modeling,
more specifically, individual-based modeling. Notably, there are various
computational modeling methods have been proposed in recent years to model
biofilm. For example, one can use molecular dynamics simulations to model the
biochemical properties of biofilm on the molecular scale Powell et al. (2018),
using dissipative particle dynamics to model biofilm deformation under shear
flow Xu et al. (2011) and coarse-grained molecular dynamics to study dewetting
phenomena Brandani et al. (2015) to model biofilm in the mesoscale and use
finite element methods to simulate the linearized growth Smith et al. (2007)
to model biofilm on the continuum scale. Among all these different methods, we
choose the individual-based modeling (IbM) Li et al. (2019) method to
represent each bacteria cell as individual particles, as it is capable of
modeling the growth and dynamics of biofilm capturing both the cellular to
clustered scales with relatively low computational resources required. More
specifically, there are three main reasons for choosing IbM:

  * •

IbM is a general multiscale method, capable of capturing the scaling effects
from cell to “film”. Since each bacteria cell is considered as a particle in
the simulation, both “cell-cell”111many cells consist the film, “cell-cell”
means the dynamics of individuals within the films, “film-materials”, and
“cell-materials” interactions are captured and described. When studying the
transport of biofilm within porous regions, where both the individual and
group dynamics play important roles, the ability to capture multiscale
mechanics is essential in our problem. The flexibility of tuning parameters
among scales is essential Li et al. (2019).

  * •

The IbM method is physically realistic to the particular scale of interest.
The scale we are focusing on is mostly the micrometer scale, where IbM offers
extremely high representation power and accuracy. First, the biofilms observed
in natural “pore-scale” mostly refer to the scale of $10^{-5}\sim 10^{-3}$m
Kapellos et al. (2015), which is particularly suitable for IbM. Note that each
bacteria cells are approximately 1$\mu$m, making the pore-scale perfectly
captures the local morphology of the biofilms. Second, the adhesion and other
micromechanisms that governs the overall mechanical behavior of biofilm mainly
originate in the micrometer scale Galy et al. (2012), where IbM could bring in
the computational tools for decent understanding. Third, our ultimate goal is
to bring in our theoretical predictions and understanding to experimental
implementations for ELM that can be used for designing new materials and
devices. Most recent work on ELM is on the micrometer scale Rodrigo-Navarro et
al. (2021), where our computational efforts could bring in the most help and
impact.

  * •

Compared with other methods, IbM has the most decent computational burden
requirement for relatively high fidelity222The goal is to combine simulation
with optimization, where the simulation is treated as the evaluated function.
Hence, the function evaluation time is important for efficient optimizations..
On the molecular scale, if one is to simulate the growth of biofilm using
methods like molecular dynamics (MD) simulation or Monte Carlo sampling, the
computational burden will be extremely high, making characterizing the
mechanism on the micrometer scale impossible. As a reference, it would require
6 months to run MD of a protein structure for 1 $\mu$s Li et al. (2021),
making this method infeasible for our problem. For the continuum scale, the
simulation of biofilm usually incorporates extended finite element methods
(XFEM) and level set method (LSM) Duddu et al. (2008), which is also extremely
computationally burdensome. To elaborate, the incorporation of the FEM with
LSM usually requires a moving mesh that resolves the phase boundary Zhai et
al. (2022), which significantly increases the computational resources
required.

To solve the second problem, our answer is to use approximation methods to
solve the inverse problem of materials design. If one is to define designing
materials by perturbing their original structures to obtain the target
properties as a forward problem, one can then define obtaining the tailored
materials’ structures from the predefined targeted properties as an inverse
problem333The rigorous formulation follows the Hadamard’s principles, which we
do not discuss in details here. The detailed inverse problem here is
formulated as finding the optimal porous structure corresponding to the target
biofilm transport properties (i.e., maximize biofilm growth), as a class 2
inverse problem. If one is to examine this defined problem in detail, there
are two main problems: (A) The defined inverse problem is ill-posed Hadamard
(1902). Two or more different porous structures may yield the same biofilm
transport properties, in which if one yields back the materials structure
representation as the solution of the inverse problem, this solution may not
be unique. (B) There are no analytical (or symbolic) forms of the inverse map.
The biofilm simulation is constituted of iterative growth and update of
bacteria cells, where it is almost impossible to obtain an analytical inverse
of this coupled multiphysics system with changed parameters444For the detail
of the simulation algorithms please refer to Refs. Li et al. (2019) & Zhai and
Yeo (2022).

To solve problem (A), our proposed solution is characterizing the design
space. We hope to approximate a surrogate model of the design space and verify
the approximated map by conducting verification simulations along the observed
maximal solution and randomly selected points. This allows us to verify the
accuracy of the fitted surrogate and the further analysis is reliable. To
solve the problem (B), our way is to avoid gradient-based optimizations and
use machine learning (ML) techniques (i.e., Gaussian process regression
(GPR)), which also allows us to do direct surrogate modeling of the design
space, and solve the first problem simultaneously. Hence, combining the two
solutions we proposed, Bayesian optimization (BO) Frazier (2018) seems to be
the perfect fit, which mainly consists of a GPR to approximate the design
space map and an acquisition function to update the solution search scheme. To
summarize and elaborate in detail, there are three major reasons for choosing
BO:

  * •

The flexibility of handling complex problems. Compared with gradient-based
methods, BO is flexible and can be adapted to solving complicated optimization
problems without requiring the calculation of the derivative of the evaluated
functions.

  * •

It is less computationally burdensome compared with other ML methods. As a
non-parametric method, GPR requires less computational resources compared with
neural networks (NN) and is especially suitable for problems defined within
the limited data regime Fuhg and Bouklas (2022). Compared with the widely used
deep reinforcement learning (DRL) Sutton and Barto (2018), BO does not require
iterative training of deep NN for each function evaluation, and hence is
significantly less computationally burdensome.

  * •

The approximation of the design space map allows direct characterization and
analysis of the sampling process. Compared with metaheuristic methods such as
genetic algorithms Mitchell (1998) or particle swarm optimization Kennedy and
Eberhart (1995), in which the updates of the function evaluation are based on
random perturbation of the input variables inspired by natural phenomena, the
learning of the design space map from GPR allows us to do detailed
characterization, and hence support our proposed solution to the problem (A).
Other than that, BO usually does not heavily rely on population, usually
requiring one evaluation per iteration, making the characterization much
easier.

In a more detailed sense, in this paper, we are combining IbM and BO to solve
a focused problem: Inversely design the porous structural materials for
biofilms transport and characterize the biomechanics from the optimization
processes. By solving this problem, we hope to potentially answer the
following questions: (1) What are the optimal porous microstructures that can
maximize the transportability of biofilms? (2) Are the approximated design
space accurate and how do we verify them? (3) What biomechanical mechanisms
are discovered by conducting the optimization and characterizing the design
space? We will answer these questions in the following sections.

This paper is formulated as follows. In Section 2 we briefly introduce the
method we used, including our computational models of biofilm physics (Section
2.1), the Bayesian optimization scheme (Section 2.2), including the surrogate
modeling of GPR (Section 2.2.1) and the iterative update scheme by acquisition
function (Section 2.2.2), followed by our designed three numerical
experimental formulations for different porous material in Section 2.3. We
then showed our results in Section 3: we discuss our optimization processes
and optimal structures for different numerical experiments in Sections 3.1 &
3.2, verify the discovered new phenomena, and provide additional mechanistic
explanations in Section 3.3. Eventually, we conclude the paper in Section 4.

##  2 Methods

As elaborated in Section 1, we will use computational methods to model the
growth of biofilms and their mechanical interactions with the porous
metamaterials in a predefined simulation box. We then combine the material
representation of the porous structure parameterized based on our defined
numerical experiments and combine the simulation framework with the Bayesian
optimization methods for iterative searching for a better porous structure
with better biofilm transport properties.

The general schematic of this study is represented in Figure 1: one begins
with an inspiration from a natural phenomenon: biofilms mostly grow in
confined systems Friedlander et al. (2013), then one can define a porous
structure that allows the biofilm to grow within to mimic this phenomenon
(Figure 1 A). What follows is one can then run the simulation initiated by
parametrized materials representation (Figure 1 B) for coupling with Bayesian
optimization (Figure 1 C). The coupling is enabled by “variable passing”
between the simulation and optimization: the simulation takes the materials’
representation as input and outputs the biofilm transport properties as the
objective to input to the optimization, where the optimization algorithm
updates the new materials’ representation as its output for an iterative loop.
This iterative searching will eventually propose an optimal structure (Figure
1 D). By characterizing the design space obtained in the optimization (Figure
1 C) and comparing the observation from these simulations, one can then
propose explanations for the optimal structures and identify new mechanisms of
biofilm transport physics (Figure 1 E).

In the following subsections, we first briefly introduce the basic formulation
of our computational methods of individual-based modeling, and then the basic
mathematical formulation of Bayesian optimization. Eventually, we briefly
introduce the designed formulation of different formulated problems for porous
membranes, lattice metamaterials, and nonconvex porous media, respectively.

Figure 1: The overall schematic of this research. The formulation is inspired
by the simulation of biofilm transport in porous materials in (A) (on top),
which inspires us to define a computational framework on the bottom, where
there are some initial bacteria cells distributed at the bottom of simulations
growing to throw a porous media as indicated in the grey area. (B) The growth
processes of the biofilms within the porous materials, i.e., the computational
simulations, calculated by individual-based modeling. (C) The reconstructed
design space for the porous material from the Bayesian optimization. (D) The
extracted optimal design from the design space. (E) One may uncover new
physical phenomena and the mechanism of bacteria transport in porous materials
by comparing the optimal design and designed benchmark cases.

###  2.1 Computational Models

In this work, our employed IbM computational models developed based on the
Newcastle University Frontiers in Engineering Biology (NUFEB) framework Li et
al. (2019), in which each bacteria cell is modeled as a spherical particle.
Biofilms are formed by cell division and extrusion of EPS. Following our
previous work on surface shape optimization Zhai and Yeo (2022), the microbe
growth and decay are governed by the following differential equation:

$\frac{dm_{i}}{dt}=\xi_{i}m_{i}$ (1)

where $m_{i}$ is the biomass of the $i^{th}$ bacteria cells and $\xi_{i}$ is
the growth rate. The growth rate of the bacteria cell is $\xi=0.00028s^{-1}$.
To avoid overlap of the particles during the growth processes, the particles
are mechanically relaxed using the individual-based approach, solved via
Newton’s equation

$m_{i}\frac{d{\bf v}_{i}}{dt}=\mathbf{F}_{c,i}+\mathbf{F}_{a,i}$ (2)

where $\mathbf{v}_{i}$ is the particles’ velocity. The contact force
$\mathbf{F}_{c,i}$ is a pair-wise force between particles to prevent
overlapping based on Hooke’s law

$\mathbf{F}_{c,i}=\sum_{j=1}^{N_{i}}\left(K_{\mathbb{N}}\delta\mathbf{n}_{i,j}-m_{i,j}\gamma_{\mathbb{N}}\mathbf{v}_{i,j}\right)$
(3)

where $N_{i}$ is the total number of neighboring particles of $i$,
$K_{\mathbb{N}}$ is the elastic constant for normal contact,
$\delta\mathbf{n}_{ij}$ is the overlap distance between the center of particle
$i$ and its neighbour particle $j$. $\gamma_{\mathbb{N}}$ is the viscoelastic
damping constant for normal contact, and $v_{i,j}$ is the relative velocity of
the two particles. The EPS adhesive force ${\bf F}_{a,i}$ is a pair-wise
interaction modelled as a van der Waals force

$\mathbf{F}_{a,i}=\sum_{j=1}^{N_{i}}\frac{H_{a}r_{i,j}}{12h_{min,i,j}^{2}}\mathbf{n}_{i,j}$
(4)

where $H_{a}$ is the Hamaker coefficient, $r_{i,j}$ is the effective outer-
radius of the $i^{\rm th}$ and $j^{\rm th}$ particles. $h_{min,i,j}$ is the
minimum separation distance of the two particles, and $\mathbf{n}_{i,j}$ is
the unit vector from particle $i$ to $j$.

Mechanical equilibrium is achieved when the average pressure of the microbial
community reaches a plateau. The average pressure $P$ of the system calculates

$P=\frac{1}{3V}\left(\sum_{i=1}^{N}m_{i}\mathbf{v}_{i}\cdot\mathbf{v}_{i}+\sum_{i=1}^{N}\sum_{j>i}^{N}\mathbf{r}_{i,j}\cdot\mathbf{F}_{i,j}\right)$
(5)

where $V$ is the sum of the particles’ volumes. The first term in the bracket
is the contribution from the kinetic energy of each particle. The second term
is the interaction energy, where $\mathbf{r}_{i,j}$ and $\mathbf{F}_{i,j}$ are
the distance and force between two interacting particles $i$ and $j$,
respectively.

Here, we employ the Monod-based method Monod (1949) to model microbial growth,
in which the growth rate is determined by the Monod kinetic equation driven by
the local concentration of nutrients. The porous materials are modeled as
fully rigid particles with neither growth nor decay. Here, under the Monod
model formulation, each bacteria cells first grow with increasing radii, and
after their radii reach a critical value $r^{\sf C}=1.36\times 10^{-6}$m, the
cell is separated into two daughter cells (Details see Ref. Li et al. (2019)).
The EPS, also modeled as particles, are secreted by the main bacteria cells in
the growing process (Details see Refs. Xavier et al. (2005); Jayathilake et
al. (2017)). After certain iterations, the system preserves a total number of
bacteria cells and the EPS particles, $\mathcal{N}_{\rm bio}^{\rm total}$.

###  2.2 Bayesian Optimization

The goal of optimization is to minimize or maximize an objective function,
which in our case is the bacteria cell number under a target design region,
denoted as $\mathcal{N}_{\rm bio}$ for ease of notation ($\mathcal{N}_{\rm
bio}\subset\mathcal{N}_{\rm bio}^{\rm total}\in\mathbb{Z}$). Using
$\mathcal{N}_{\rm bio}=\mathcal{M}_{\sf NUFEB}(N_{\rm
unit},\bar{\mathcal{D}};\mathbf{p})$ to denote a multivariate function
relation, in which $N_{\rm unit}$ and $\bar{\mathcal{D}}$ stand for unit cell
numbers per simulation box side and the dimensionless structural parameter (or
dimensionless variable), respectively, are the design variables to be
elaborated in details in Section 2.3. For simplicity, we use
$\mathcal{DV}=[N_{\rm unit}$, $\bar{\mathcal{D}}]$ to denote the design
variables. $\bf p$ is the parameter involved in the numerical simulation, as
presented in Equations (1$\sim$5). The optimization process can be simplified
as:

$\displaystyle\operatorname*{arg\,max}_{N_{\rm
unit},\bar{\mathcal{D}}}\mathcal{N}_{\rm bio}={\mathcal{M}_{\sf NUFEB}}(N_{\rm
unit},\bar{\mathcal{D}};\mathbf{p}),$ (6) $\displaystyle{\rm subject\
to}\quad{\bar{\mathcal{D}}}_{\rm
LB}\leq\bar{\mathcal{D}}\leq\bar{\mathcal{D}}_{\rm UB},\ 1\leq N_{\rm
unit}\leq 15\ (N_{\rm unit}\in\mathbb{Z})$

Here, we define a target growth region to count $\mathcal{N}_{\rm bio}$
(Section 2.3), so that the optimizations are tailor the materials’
microstructure to enhance growth toward the target region. Given the input
design variables $\mathcal{DV}$, we represent the biofilm physics growth
simulation model as a map, $\mathcal{M}_{\sf NUFEB}:N_{\rm
unit},\bar{\mathcal{D}}\rightarrow\mathcal{N}_{\rm bio}$, where the simulation
parameters
$\mathbf{p}=[\xi_{i},K_{\mathbb{N}},\gamma_{\mathbb{N}},H_{a},r^{\mathsf{C}},...]$
are incorporated in the IbM model (Section 2.1). $\mathcal{M}_{\sf
NUFEB}(\cdot)$ stands for the numerical simulation from NUFEB that maps the
design representation of the materials as input and the bacterial cell number
count as output. $N_{\rm unit}$ is an integer between 1 and 15 as the number
of unit cells are changing along the BO iterations. The dimensionless
structure parameter $\bar{\mathcal{D}}$ is defined per case, as the lower and
upper bounds $\bar{\mathcal{D}}_{\rm LB}$ & $\bar{\mathcal{D}}_{\rm UB}$
differs based on the simulation and materials basis settings, to be discussed
in Section 2.3.

BO aims to iteratively update new evaluations from the computational models in
Section to search for optimal materials. Through sampling multiple simulations
and mapping the design variables into the defined objective, one can construct
a surrogate of the direct map between the input (i.e., the design variables)
and the output (i.e., the objective) from GPR. This GPR reconstructed
surrogate is then updated through the acquisition functions of choice.

####  2.2.1 Gaussian Process Regression

GPR is a Bayesian statistical approach to approximate and model function(s).
Considering our optimization problem, the function can be denoted as
$\mathcal{N}_{\rm bio}=\mathcal{M}_{\sf NUFEB}(\mathcal{DV};\mathbf{p})$,
where $\mathcal{N}_{\rm bio}$ is evaluated at a collection of different sets
of points (or design variables):
$\mathcal{DV}_{1},\mathcal{DV}_{2},...,\mathcal{DV}_{k}\in\mathbb{R}^{2}$, we
can obtain the vector $[\mathcal{M}_{\sf
NUFEB}(\mathcal{DV}_{1}),...,\mathcal{M}_{\sf NUFEB}(\mathcal{DV}_{k})]$ to
construct a surrogate model for the design parameters with the correlated
objectives. The vector is randomly drawn from a prior probability
distribution, where GPR takes this prior distribution to be a multivariate
normal with a particular mean vector and covariance matrix. Here, the mean
vector and covariance matrix are constructed by evaluating the mean function
$\mu_{0}$ and the covariance function $\Sigma_{0}$ at each pair of points
$\mathcal{DV}_{i}$, $\mathcal{DV}_{j}$. The resulting prior distribution on
the vector $[\mathcal{M}_{\sf NUFEB}(x_{1}),...,\mathcal{M}_{\sf
NUFEB}(x_{k})]$ is represented in the form of a normal distribution to
construct the surrogate model Frazier (2018):

$\mathcal{N}_{\rm
bio}(\mathcal{DV}_{1:k})\sim{\mathfrak{N}}\left(\mu_{0}(\mathcal{DV}_{1:k}),\Sigma_{0}(\mathcal{DV}_{1:k},\mathcal{DV}_{1:k}))\right)$
(7)

where $\mathfrak{N}(\cdot)$ denotes the normal distribution. The collection of
input points is represented in compact notation: $1:k$ represents the range of
$1,2,...,k$. The surrogate model $\mathcal{M}_{\sf NUFEB}(\mathcal{DV})$ on
$1:k$ is represented as a probability distribution given in Equation (7). To
update the model with new observations, such as after inferring the value of
$\mathcal{M}_{\sf NUFEB}(\mathcal{DV})$ at a new point $\mathcal{DV}$, we let
$k=l+1$ and $\mathcal{DV}_{k}=\mathcal{DV}$. The conditional distribution of
$\mathcal{N}_{\rm bio}$ given observations $\mathcal{DV}_{1:l}$ using Bayes’
rule is

$\displaystyle\mathcal{N}_{\rm bio}(\mathcal{DV})|\mathcal{N}_{\rm
bio}(\mathcal{DV}_{1:l})$
$\displaystyle\sim\mathfrak{N}(\mu_{l}(\mathcal{DV}),\sigma_{l}^{2}(\mathcal{DV}))$
(8) $\displaystyle\mu_{l}(\mathcal{DV})$
$\displaystyle=\Sigma_{0}(\mathcal{DV},\mathcal{DV}_{1:l})\Sigma_{0}(\mathcal{DV}_{1:l},\mathcal{DV}_{1:l})^{-1}\left(\mathcal{M}_{\sf
NUFEB}(\mathcal{DV}_{1:l})-\mu_{0}(\mathcal{DV}_{1:l})+\mu_{0}(\mathcal{DV})\right)$
$\displaystyle\sigma_{l}^{2}$
$\displaystyle=\Sigma_{0}(\mathcal{DV},\mathcal{DV})-\Sigma_{0}(\mathcal{DV},\mathcal{DV}_{1:l})\Sigma_{0}(\mathcal{DV}_{1:l},\mathcal{DV}_{1:l})^{-1}\Sigma_{0}(\mathcal{DV}_{1:l},\mathcal{DV})$

where the posterior mean $\mu_{l}(\mathcal{DV})$ is a weighted average between
the prior $\mu_{0}(\mathcal{DV})$ and the estimation from $\mathcal{M}_{\sf
NUFEB}(\mathcal{DV}_{1:l})$, where the weight applied depends on the kernel
used.

Here, we use the Gaussian kernel, hence the prior covariance is Biswas et al.
(2021)

$\displaystyle\Sigma_{0}(\mathcal{DV}_{i},\mathcal{DV}_{j})$
$\displaystyle=\sigma^{2}R(\mathcal{DV}_{i},\mathcal{DV}_{j}),$ (9)
$\displaystyle R(\mathcal{DV}_{i},\mathcal{DV}_{j})$
$\displaystyle=\mathtt{exp}\left(\frac{1}{2}\sum_{m=1}^{d}\frac{\left(\mathcal{DV}_{i,m}-\mathcal{DV}_{j,m}\right)^{2}}{\theta_{m}^{2}}\right)$
$\displaystyle\theta_{m}$
$\displaystyle=(\theta_{1},\theta_{2},...,\theta_{d})$

where $\sigma^{2}$ is the overall variance parameter and $\theta_{m}$ is the
correlation length scale parameter in dimension $m$ of the $d^{\rm th}$
dimension of $\mathcal{DV}$, which are all hyperparameters of GPR.
$R(\mathcal{DV}_{i},\mathcal{DV}_{j})$ is the spatial correlation function.
Our goal is to estimate the parameters $\sigma$ and $\theta_{m}$ that create
the surrogate model given the training data $[(\mathcal{N}_{\rm bio})_{k},\
\mathcal{DV}_{k}]$ at iteration $k$. Here, we will use $\hat{\mathcal{M}}_{\sf
GPR}$ to denote the surrogate model constructed from GPR in the iterative
updating process. The updating sampling scheme is achieved through the
acquisition function in the following section, which improves the accuracy of
the updated surrogate so that the reconstructed design space is approximating
the theoretical continuous design from NUFEB simulations
$\hat{\mathcal{M}}_{\sf GPR}\sim\mathcal{M}_{\sf NUFEB}$.

####  2.2.2 Acquisition Function

Given the training data $[(\mathcal{N}_{\rm bio})_{k},\ \mathcal{DV}_{k}]$,
Equation (7) gives us the prior distribution $(\mathcal{N}_{\rm
bio})_{l}\sim\mathfrak{N}(\mu_{0},\Sigma_{0})$ as the surrogate. This prior
and the given dataset induce a posterior: the acquisition function denoted as
$\mathcal{A}:\mathcal{X}\xrightarrow{}\mathbb{R}^{+}$, determines the point in
$\mathcal{X}$ to be evaluated through the proxy optimization
$\mathcal{DV}_{\rm
best}=\operatorname*{arg\,max}_{\mathcal{DV}}\mathcal{A}(\mathcal{DV})$. The
acquisition function depends on the previous observations, which can be
represented as
$\mathcal{A}=\mathcal{A}(\mathcal{DV};(\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}),\theta)$, where $(\mathcal{DV}_{l},(\mathcal{N}_{\rm bio})_{l})$
leads to the reconstructed $\hat{\mathcal{M}}_{\sf GPR}$. Taking our previous
notation, the new observation is probed through the acquisition Deshwal et al.
(2021):

$\mathcal{DV}_{k}=\mathcal{DV}_{l+1}=\operatorname*{arg\,max}_{\mathcal{DV}}\mathcal{A}\left(\mathcal{DV};(\hat{\mathcal{M}}_{\sf
GPR})_{l},\theta_{m}\right)$ (10)

where the input space contains the evaluation of design variables at $l$
points: $(\mathcal{DV}_{1},\mathcal{DV}_{2},...,\mathcal{DV}_{l})$. We compare
and characterize two different acquisition functions, the Upper Confidence
Bound (UCB) and the Expected Improvement (EI), to do a benchmark study on the
effect of acquisition updates. The UCB exploits the upper confidence bounds to
construct the acquisition and minimize the regret. UCB takes the form Snoek et
al. (2012)

$\mathcal{A}_{\sf UCB}\left(\mathcal{DV};(\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}),\theta_{m}\right):=\mu_{l}\left(\mathcal{DV};(\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}),\theta_{m}\right)+\kappa\sigma\left(\mathcal{DV};(\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}),\theta_{m}\right)$ (11)

where $\kappa$ is a tunable parameter balancing exploitation and exploration
when constructing the surrogate model. We take $\kappa=2$ in our
implementations. For the EI acquisition, the function writes:

$\displaystyle\mathcal{A}_{\sf
EI}\left(\mathcal{DV};(\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}),\theta_{m}\right):=\sigma_{l}\left(\mathcal{DV};(\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}),\theta_{m}\right)\left(\gamma(\mathcal{DV}){\bf\Phi}\left(\gamma(\mathcal{DV})\right)+\mathfrak{N}\left(\gamma(\mathcal{DV});0,1\right)\right)$
(12)

where $\gamma$ is computed as $\gamma=\left(-\mathcal{M}_{\sf
NUFEB}(\mathcal{DV}_{\rm
best})+\mu(\mathcal{DV};\\{\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}\\}_{l},\theta)-\Xi\right)/\sigma\left(\mathcal{DV};\\{\mathcal{DV}_{l},(\mathcal{N}_{\rm
bio})_{l}\\}_{l},\theta\right)$, where $\Xi$ is a damping factor in the code
implementation, and $\Xi=10^{-4}$ in our implementation. Note that
$\mathcal{A}_{\sf EI}$ preserves a closed form under the GP evaluations.

Combining GPR and the acquisition function, the surrogate model can
approximate the design space’s maximal value. In our case, such BO methods are
applied to obtain optimal porous materials’ structures with maximal bacterial
cell numbers in the target transport region. Here, the total function
evaluations are different per case, as to be discussed in the following
Section 2.3.

###  2.3 Numerical Experiments

Figure 2: The schematic illustration for the three different porous materials
formulations. The porous materials are treated as repeated elements of unit
cells, and the number of unit cells per length is $N_{\rm unit}$ (marked in
the middle sub-figures), which is defined as a design variable in the
optimization. For every unit cell, the dimensionless structure parameter,
$\bar{\mathcal{D}}$, is defined to quantify the vacuum-solid region spatial
ratio in a defined unit cell illustrated in the right sub-figures. (A) Two-
dimensional porous membranes for biofilm transport. Note that the “two-
dimensional” means no repeated unit cells are in the third dimension, i.e.,
the Z axis. The design variables hence do not perturb the geometries in the
third dimension. Bacteria cells are grown within the “micro-pipelines” within
the membranes to the top region. The dimensionless variable writes
$\bar{\mathcal{D}}=\mathcal{R}_{\rm vol}/(\mathcal{R}_{\rm
vac}+\mathcal{R}_{\rm vol})$, is defined as the radii ratio between the vacuum
region and the overall region (vacuum + volumetric solid). (B) Lattice porous
metamaterials for biofilm transport. Bacteria cells are grown within the
porous region within the lattice microstructures to reach the top. The unit
cell dimensionless variable takes the form $\bar{\mathcal{D}}=\ell_{\rm
vac}/(\ell_{\rm vac}+\ell_{\rm vol})$, is defined as the length ratio between
the vacuum region and the overall region. (C) Non-convex three-dimensional
porous media for biofilm transport. Bacteria cells are grown within the porous
region within the porous media to reach the top. The unit cell dimensionless
variable takes the form $\bar{\mathcal{D}}=\mathcal{R}_{\rm
vol}/(\mathcal{R}_{\rm vac}+\mathcal{R}_{\rm vol})$, is defined as the radii
ratio between the volumetric region and the overall region.

Here, we define three different simulation cases to simulate the process of
biofilm growth constrained within porous materials, inspired by experimental
setup, literature results, and natural phenomena. The general schematic
representing the numerical experimental setup is illustrated in Figure 2.
Recall the overall optimization formulation in Equation (6), $\mathcal{N}_{\rm
bio}$ is the bacteria cell numbers counted in the top quarter region nominated
as the objective growth region, i.e., $L_{\rm obj}\times L_{X}\times L_{Y}$.
The porous materials’ microstructures are defined in the materials region,
i.e., $L_{\rm mat}\times L_{X}\times L_{Y}$. The initial bacteria cells are
distributed in the initial biomass region, i.e., $L_{\rm bio}\times
L_{X}\times L_{Y}$. $N_{\rm unit}$ are formulated differently based on the
“dimension” of the problem, where for the porous membrane (Figure 2 A) $N_{\rm
unit}$ is only defined in the X-Y plane. For lattice metamaterials and non-
convex porous media, it is defined in both the X, Y, and Z directions.
$\bar{\mathcal{D}}$ are defined within the unit cells. Here,
$L_{X}=L_{Y}=50\mu\rm m$, $L_{\rm bio}=L_{\rm obj}=12.5\mu\rm m$, and $L_{\rm
mat}=25\mu\rm m$.

  * •

Porous Membranes. Considering biofilm growth and flow constrained in a
microchannel are widely applied and studied by the microfluidics communities
and their wide applications spanning from energy, biosensing, and many others
Pousti et al. (2019); Ye et al. (2021); many numerical Landa-Marbán et al.
(2019); Aspa et al. (2011) and theoretical Landa-Marbán et al. (2020) studies
have been developed trying to understand the mechanism of biofilm growth and
flow in a microchannel. Here, our numerical implementations for channeled
biofilm growth are mainly inspired by the simulation setup by Aspa et al. Aspa
et al. (2011), where cylinder-shaped convex pores are “drilled” in the solid
materials to create channels for biofilm to grow within in Figure 2 A. The
morphology of the unit cell is shown in the right subfigure in Figure 2 A: the
radius of the hole (vacuum area) is denoted as $\mathcal{R}_{\rm vac}$ and the
length of the residual solid body (the volumetric part, equals to half length
of the unit cell minus $\mathcal{R}_{\rm vac}$) is denoted as
$\mathcal{R}_{\rm vol}$. The dimensionless variable can then be computed as
$\bar{\mathcal{D}}=\frac{\mathcal{R}_{\rm vac}}{\mathcal{R}_{\rm
vac}+\mathcal{R}_{\rm vol}}$. In this scenario, the range of the dimensionless
variable is defined as $\bar{\mathcal{D}}\in[0.1,0.9]$
($\bar{\mathcal{D}}_{\rm LB}$ and $\bar{\mathcal{D}}_{\rm UB}$ in Equation
(6)). The meaning of defining this problem is the optimization results from
designing porous channels (or 2D porous membranes) could potentially bring in
solutions for biofilm transport and utilization as ELM, as this kind of
topological formulation is easy to manufacture. Also, based on this material
formulation, we also conduct a benchmarking study comparing the effect of
acquisition function in sampling the design space from BO (Section 2.2.2), in
which we also characterize the design space from the sampling perspective,
that could guide general materials design optimizations.

  * •

Lattice Metamaterials. In recent five years, there is a huge growth in the
study of the design Ma et al. (2022); Shaw et al. (2019) and properties Gu
(2018); Portela et al. (2020) of mechanical metamaterials recently (or
synonymously architectured materials). However, their potential applications
in biomass storage and transport are rarely touched, with very few works
concerning their potential use as biofilm carriers Ovelheiro (2020); He et al.
(2021) and related properties Hall et al. (2021). Here, we hope to use our
simulations to fill in this gap and brought new insights into the
possibilities of using lattice metamaterials for biofilm storage and
transport. The unit cell of such metamaterials is shown in the right subfigure
Figure 2 B: the half length of the vacuum area is denoted as $\ell_{\rm vac}$,
and the edge length of the solid volumetric part is denoted as $\ell_{\rm
vol}$, where the dimensionless variable is defined as
$\bar{\mathcal{D}}=\frac{\ell_{\rm vol}}{\ell_{\rm vac}+\ell_{\rm vol}}$. The
range of the dimensionless variable is defined as $[0.1,0.5]$.

  * •

Non-convex Porous Media. Inspired by the fact that biofilms were mostly found
in natural habitats where they were constrained in pseudo- or spherical solid
bodies Bhattacharjee and Datta (2019); Carrel et al. (2018); Coyte et al.
(2016); Kurz et al. (2022), we propose the simulation scenario where biofilm
grows in nonconvex solid bodies shown in Figure 2 C. The simulations were
mainly inspired by the study of Dehkharghani et al. Dehkharghani et al. (2023)
and Bhattacharjee & Datta Bhattacharjee and Datta (2019), where we are using
BO as a tool to sample the scale effect studied in Dehkharghani et al. (2023)
defined a similar 3D porous packing of solid spherical bodies in Bhattacharjee
and Datta (2019). The dimensionless variable is defined as the radii ratio
between the solid spheres and the overall unit cell lengths (right subfigure
in Figure 2 C): $\bar{\mathcal{D}}=\frac{\mathcal{R}_{\rm
vol}}{\mathcal{R}_{\rm vol}+\mathcal{R}_{\rm vac}}$. The range of the
dimensionless variable is defined as $[0.5,1.2]$. Note that instead of
simulating potential manufacturable porous materials to inspire industrial
applications, we hope to use this case to combine with the BO sampling to
investigate the biofilm transport scenarios in nature.

We use the porous membrane case to first characterize the acquisition function
used, and apply the BO for 500 iterations each. For the lattice metamaterials
case, due to the high computation burden of the simulation, we only apply BO
for 300 iterations with only the EI acquisition function. For the porous media
case, we apply BO for 500 iterations with only the EI acquisition function.
For both three cases, we conduct characterization simulation cases to examine
the accuracy of the GPR approximated design space, in which one case is based
on observation towards the maximal point in the visualized reconstructed
design space, and the other case is based on random test selected in the
design space.

##  3 Results & Discussion

###  3.1 Porous Membranes

Figure 3 shows the whole optimization process updated by both EI and UCB
acquisition functions for the porous membrane design case. Figure 3 A
visualizes the change of the objective along with the iteration, in which the
red line stands for the optimization process updated by the EI acquisition
function, and the blue line stands for that of the UCB acquisition function.
It can be observed that UCB exhibits more evident fluctuation along the
sampling process and the EI acquisition sampled objectives are more
“clustered” towards the upper region. To be more rigorous, we generate Figure
3 B, visualizing the overall statistical distribution of the objectives by two
different acquisition functions. It can be qualitatively observed that the
variance of EI is evidently smaller than that of UCB, and the mean objective
value sampled by EI is higher. Quantitatively, the objective variance for the
EI and UCB acquisition functions are $2.62\times 10^{7}$ and $3.46\times
10^{7}$, respectively, where the UCB acquisition sampled objectives’ variance
is relatively 32.08% higher than the EI acquisition. The mean objective values
updated by EI and UCB acquisition functions are $\mathcal{N}_{\rm bio}^{\sf
EI}=30502$ & $\mathcal{N}_{\rm bio}^{\sf UCB}=30056$, respectively. The EI
mean objective is relatively 1.49% higher than the UCB acquisition function.
Figure 3 C1 & C2 visualizes the trends of the normalized design variables
along the sampling process by EI and UCB acquisition functions, respectively.
It can be deduced from both the subfigures that $\bar{\mathcal{D}}$ are
generally being sampled towards higher values and $N_{\rm unit}$ are being
sampled in relatively lower values during the optimization processes, by
observing their value range visualized by the color bar.

Figure 3: Optimization results for 2D porous membrane design. (A) The change
of the objective value $\mathcal{N}_{\rm bio}$ along the iteration process.
The red dotted dashed line represents the BO process through the Expected
Improvement (EI) acquisition function. The blue dotted dashed line represents
the BO process through the Upper Confidence Bound (UCB) acquisition function.
(B) The statistical distribution of the objective values along the
optimization processes characterized by the two different acquisition
functions. (C) The normalized design variable change along the iteration
process, corresponding to subfigure A, where subfigures C1 & C2 represent the
BO updated by EI and UCB acquisition function, respectively.

Based on our qualitative observations from Figure 3, three questions may
naturally arise as further verifying the qualitative observations.

  1. (I)

Just observing the objectives change may not be comprehensive enough to
estimate whether both the acquisition functions are sampling toward the
“correct” directions, i.e., whether the sampling directions are moving toward
higher objective values (i.e., the design goal).

  2. (II)

Can we generally verify the accuracy of the approximated GPR approximated
design space?

  3. (III)

Concerning the observations in Figure 3 C1 & C2, what are the exact geometries
represented by the changing variables?

Note that these three questions are fundamental in our following analysis for
different materials design cases. Here, to answer Question (I), we generate
Figure 4, to visualize the sampling process during the optimizations, and
characterize them with the overall sampling density. To answer Questions (II)
& (III), we generate Figure 5, to characterize the approximated design space
using simulations and visually show the general trends captured by the
approximated models and simulation points. We then further visualize the
geometries extracted from the characterization simulations.

Figure 4: The design space reconstruction (visualized in normalized values)
and sampling density maps by the two different acquisition functions for 2D
porous membrane design case. Here, (A1$\sim$C1) stand for the design space
surrogate and sampling density map from the EI acquisition function, and
(A2$\sim$C2) stand for those by UCB acquisition function. Note that for
subfigures A, the white dots are visualized in three batches: the first batch
represents the first 300 iterations, visualized in small circular dots, the
mid-100 iterations are visualized as squared-shaped dots, and the last 100
iterations are visualized in large triangular dots, which are the easiest to
be identified. For subfigure B, the visualization of the first two batches
remains the same, whereas the last batch set contains different evaluations
and is marked still in triangular dots. For details please see the texts. The
main goal is to characterize the sampling density map through the morphology
of the sampling dots in the reconstructed design space. (A1) The reconstructed
design space by EI acquisition function. (B1) Zoomed view toward the target
design region from subfigure A1, where $N_{\rm unit}\in[5,15]$, and
$\bar{\mathcal{D}}\times 100\in[50,100]$. (C1) The normalized sampling density
map for the EI acquisition function, visualizing the density of the choice of
the design variables in the optimization processes. (A2) The reconstructed
design space by UCB acquisition function. (B2) Zoomed view toward the target
design region from subfigure A2, where $N_{\rm unit}\in[5,15]$, and
$\bar{\mathcal{D}}\times 100\in[50,100]$. (C2) The normalized sampling density
map for the UCB acquisition function, visualizing the design variables’
densities in the optimization processes.

Figure 4 A1 & A2 visualize the overall reconstructed design spaces updated by
EI and UCB acquisition functions. Note that the dimensionless variable
$\bar{\mathcal{D}}$ is multiplied by 100 in the visualizations for ease of
analysis. It can be observed two different acquisition functions all
approximated the same trend: there is a large objective gradient changing from
the bottom-right corner. Physically, this would indicate that when the pores’
radii ($\mathcal{R}_{\rm vac}$ in Figure 2 A) are small and the unit cell
numbers ($N_{\rm unit}$) are generally larger, the biofilm transport
capability of the porous materials decreases. One also observes that the
objective values are qualitatively higher with higher $\bar{\mathcal{D}}$
values, i.e., $\bar{\mathcal{D}}\gtrapprox 0.5$. We hence visualize the “upper
design space” in Figure 4 B1 & B2, in which the region
$\bar{\mathcal{D}}\in[0.5,0.9]$ are visualized. It can be observed that the
objective values are higher in the “top-right” corner of the design space,
where both the sampling points’ density and normalized objective values are
higher. To directly visualize the (normalized) sampling density, Figure 4 C is
created. We observe that the sampling density distribution basically overlaps
with our observations on the design space: there are higher sampling densities
toward the top-right corners (i.e., higher $\bar{\mathcal{D}}$ and $N_{\rm
unit}$ values) characterized by both acquisition functions. Combining both
Figure 4 A, B, & C, one deduces that both the reconstructed design spaces and
the sampling densities both tell us when the porous membranes contain larger
$\bar{\mathcal{D}}$ & $N_{\rm unit}$ values the biofilm transportability,
i.e., $\mathcal{N}_{\rm bio}$, are higher. Here, the EI acquisition function
samples 407 points in the “upper design space” (Figure 4 B1), and the UCB
acquisition function samples 373 points (Figure 4 B2). If we define design
space in Figure 4 B as the target region, the EI acquisition sampling
technique is 9.12% relatively more efficient than the acquisition function. If
we only look at the last 100 iterations from the BO, the EI acquisition
function samples 87 points in the target region, and the UCB acquisition
function samples 85 points. Compared with a uniformly distributed grid search
method, the EI acquisition function is 74% more efficient and the UCB
acquisition function is 70% more efficient, compared with the grid search
sampling. The EI acquisition function is 2.35% more efficient than the UCB
acquisition by estimating the last 100 design space samples in the target
region. To cross-verify these cross-validated observations from a more
quantitative perspective and answer our Questions (II) & (III), we conduct
design space characterization from additional simulations in Figure 5.

Figure 5 A & B show the general and zoomed views of the design space
characterizations comparing the selected characterization simulations (in
colored dots) and randomly selected simulations (in grey dots) to verify the
effect of the design variables ($\bar{\mathcal{D}}$ & $N_{\rm unit}$) to the
target bacteria cell numbers $\mathcal{N}_{\rm bio}$. Here, the blue dots and
grey dots in Figure 5 A are extracted based on $\bar{\mathcal{D}}=0.9$ and
0.2, respectively. The blue dots and grey dots in Figure 5 B are extracted
based on $N_{\rm unit}=15$ and 10, respectively. The $\bar{\mathcal{D}}$ and
$N_{\rm unit}$ values for blue and red dots are selected based on observations
from Figure 4 as our guess for the porous materials’ geometries that contain
the highest objective value. The $\bar{\mathcal{D}}$ and $N_{\rm unit}$ values
for the grey dots are randomly selected to compare with our observational
guess. We then directly visualize the points from the characterization
simulations on the GPR reconstructed design space in Figure 5 D: The black
triangular dots indicate the blue and red dots and the grey dots indicate the
grey dots in the left subfigures (Figure 5 A & B). It can be observed that the
characterization simulation tests fit well with the GPR-approximated design
space as both the black and grey dots overlap well with the surface contours.
We then pick a series of representative points from the characterization
simulations and directly visualize them in Figure 5 C, marked in red triangles
in Figure 5 A & B, nominated as $\mathbb{T}_{\alpha}\sim\mathbb{T}_{\gamma}$ &
$\mathbb{T}_{\rm a}\sim\mathbb{T}_{\rm c}$. It can be detected from the zoom
view in Figure 5 A that $\mathbb{T}_{\alpha}$ is evidently smaller than that
of $\mathbb{T}_{\beta}$ and $\mathbb{T}_{\gamma}$, and we can further deduce
that the porous membrane with larger pores does not necessarily enhance the
transportability of the porous materials, which is not intuitive. We propose
that the reactive forces in the pore wall drive the new bacteria cells to be
generated toward the upper region. When the radii of the pores are too large,
such reactive force does not act on the bacteria cells as they do under
smaller radii. Moreover, it can be observed from Figure 5 B that for $N_{\rm
unit}=10$ & 15, the effects of the dimensionless variable $\bar{\mathcal{D}}$
on the objective $\mathcal{N}_{\rm bio}$ are similar, where there are sudden
increases of the objective between $\bar{\mathcal{D}}\in[0.2,0.4]$.

Based on our analysis, we deduce that the EI acquisition function outperforms
the UCB acquisition function in our formulated porous materials design case by
estimating the objective variance, the mean objective values, and sampling
improvements over the design space. We also observe that with larger relative
radii of the pores and more unit cells per side, the transportability of
porous materials to biofilms is then higher, from analyzing the design space.
We then will adopt only the EI acquisition function and conduct further
analysis for lattice and 3D nonconvex porous media (Section 3.2).

Figure 5: Design space characterization for the Gaussian process regression
(GPR) reconstructed design space and topologies extraction from the
characterization processes for the 2D porous membrane design case. (A)
Characterization of the design variable $N_{\rm unit}$ with different fixed
values of $\bar{\mathcal{D}}$. Note that the blue circular dots correspond to
the black triangular dots, and the grey circular dots corresponds to the black
triangular dots, in subfigure D. The blue and red circular dots are the
characterization tests informed by qualitative observation of the GPR
reconstructed design space to approximate the optimal design (i.e., maximal
point), and the grey dots are random tests to benchmark our characterization
informed by the observations. The zoomed view describes the detailed
differences between the two sets of characterization simulations, in which
three sets of membrane topologies are selected and highlighted in red
triangular plots, nominated as $\mathbb{T}_{\alpha}$, $\mathbb{T}_{\beta}$,
and $\mathbb{T}_{\gamma}$, respectively. (B) Design variable characterization
for $\bar{\mathcal{D}}$ compared with random benchmark test marked in red and
grey dots, respectively. The zoomed view describes the detailed differences
between the two sets of characterization simulations, in which three sets of
membrane topologies are selected and highlighted in red triangular plots,
nominated as $\mathbb{T}_{\rm a}$, $\mathbb{T}_{\rm b}$, and $\mathbb{T}_{\rm
c}$, respectively. (C) Extracted porous membranes’ topologies
($\mathbb{T}_{\alpha}\sim\mathbb{T}_{\gamma}$ & $\mathbb{T}_{\rm
a}\sim\mathbb{T}_{\rm c}$) from characterizing both the design variables
$N_{\rm unit}$ and $\bar{\mathcal{D}}$ corresponding to the selections in
subfigures A & B. (D) The characterization data match with the GPR
reconstructed design spaces from both the EI and UCB acquisition function. The
black triangular dots are the characterization informed by observation from
the GPR reconstructed design space towards the maximal value. The grey
triangular dots are randomly selected test points to benchmark the
observation-informed characterizations. For details please see the texts.

###  3.2 Lattice and Porous Materials

Figure 6 shows the reconstructed design space and the sampling process along
with the sampling density updated by the EI acquisition function, similar to
what we did in Figure 4 to answer our Question (I). It can be observed from
Figure 6 A1 that the reconstructed design space from 300 evaluations is much
more nonconvex compared with that of the 2D porous membrane (Figure 4 A) and
porous media (Figure 6 A2), but the sampling are more concentrated toward the
mid-top region ($N_{\rm unit}\approx 0.5$ & $\bar{\mathcal{D}}\in[0.4,0.5]$).
Figure 6 B1 is created to better visualize this region ($N_{\rm
unit}\in[1,10]$ & $\bar{\mathcal{D}}\in[0.3,0.5]$), in which by qualitative
estimation one deduces that there are more sampling points around $N_{\rm
unit}=6$ and $\bar{\mathcal{D}}=0.5$. Comparing the reconstructed design space
and the sampling density (Figure 6 C1), one observes that the general trends
of the sampling density and the reconstructed design space overlap well, where
we thence pick $N_{unit}=6$ and $\bar{\mathcal{D}}=0.45$ for further
characterization simulations based on qualitative observations (Figure 7 1).
Figure 6 A2 shows that the reconstructed design space shapes like a “tilted
wave” — the higher objective values are distributed along the “cross-split”
across the design space coordinates. By observing both the Figure 6 A2 & C2 we
deduce that the sampling density is more centered toward the “upper design
space”. Hence, we only extract the zoomed view of the top-mid design space in
Figure 6 B2 ($N_{\rm unit}\in[1,10]$ & $\bar{\mathcal{D}}\in[0.9,1.2]$). From
Figure 6 B2 we pick $N_{\rm unit}=7$ and $\bar{\mathcal{D}}=1.1$ to conduct
characterization tests in Figure 7 2.

To estimate the effect of the acquisition function over the sampling of the
design space, we also estimate the spatial distribution of the last 100
iterations within the target design space (or target region), where the target
regions are defined based on the zoomed design space in Figure 6 B ($N_{\rm
unit}\in[1,10]\ \&\ \bar{\mathcal{D}}\in[0.3,0.5]$ for lattice metamaterials
in Figure 6 B1; $N_{\rm unit}\in[1,10]\ \&\ \bar{\mathcal{D}}\in[0.9,1.2]$ for
3D porous media in Figure 6 B2). For the lattice metamaterials, there are 62
points sampled in the target region, which is 92.89% higher than the uniform
distribution of 100 points with assumed grid search methods (32.14 points in
the target region). For the 3D porous media, there are 89 points sampled in
the target region, which is 223.04% more efficient than the uniformed sampled
100 points (27.55 points in the target region). In fine, BO exhibits
outstanding ability for sampling towards the target design goal for both
porous structures cases.

Figure 6: The design space reconstruction (visualized in normalized values)
and sampling density maps by the two different acquisition functions for
lattice metamaterials (A1$\sim$C1) and 3D porous media (A2$\sim$C2), updated
by the EI acquisition function. The morphologies of the white dots are
separated into three different batches. (A1) The reconstructed design space by
EI acquisition function. The first batch represents the first 100 iterations,
visualized in small circular dots, the mid-100 iterations are visualized as
squared-shaped dots, and the last 100 iterations are visualized in large
triangular dots. (B1) Zoomed view toward the target design region from
subfigure A1, where $N_{\rm unit}\in[1,10]$, and $\bar{\mathcal{D}}\times
100\in[30,50]$. The first batch represents the first 100 iterations,
visualized in small circular dots, the mid-50 iterations are visualized as
squared-shaped dots, and the large triangular dots represent the rest
visualizations. For details please see the texts. (C1) The normalized sampling
density map for the EI acquisition function for the lattice metamaterials
design case, visualizing the density of the choice of the design variables in
the optimization processes. (A2) The reconstructed design space by EI
acquisition function. The first batch represents the first 300 iterations,
visualized in small circular dots, the mid-100 iterations are visualized as
squared-shaped dots, and the last 100 iterations are visualized in large
triangular dots. (B2) Zoomed view toward the target design region from
subfigure A2, where $N_{\rm unit}\in[1,10]$, and $\bar{\mathcal{D}}\times
100\in[90,120]$. The first batch represents the first 300 iterations,
visualized in small circular dots, the mid-50 iterations are visualized as
squared-shaped dots, and the large triangular dots represent the rest
visualizations. For details please see the texts. (C2) The normalized sampling
density map for the EI acquisition function for the 3D porous media design
case, visualizing the design variables’ densities in the optimization
processes.

Other than our selected characterization tests, same with our porous membrane
design space characterization (Figure 5), we also randomly pick two additional
characterization tests for characterizing $N_{\rm unit}$ and
$\bar{\mathcal{D}}$ for each porous materials design case, respectively. For
designing the lattice metamaterials, we pick $\bar{\mathcal{D}}=0.1$ (Figure 7
A1) and $N_{\rm unit}=15$ (Figure 7 B1), and for 3D porous media design, we
pick $\bar{\mathcal{D}}=0.5$ (Figure 7 A2) and $N_{\rm unit}=15$ (Figure 7
B2). It can be observed from Figure 7 A & B that the selected characterization
tests generally capture the geometries of the highest objectives, where the
blue and red dots exhibit higher values than the grey dots. Interestingly, for
both porous materials cases, the topology corresponds to the highest objective
value selected from the characterization tests for $\bar{\mathcal{D}}$ (Figure
7 B), $\mathbb{T}_{\beta}$ are not the topology that contends the highest
objective value by characterizing $N_{\rm unit}$ (Figure 7 A). This indicates
that our observational guess toward the highest objective is not fully
accurate, where our characterization tests correct our initial guess and
contends the porous structural topologies $\mathbb{T}_{\alpha}$. By observing
Figure 7 D1 & D2 we observe that the characterization tests generally match
well with the GPR approximated design space, indicating the effectiveness of
the general data-driven design scheme. Notwithstanding, by comparing Figure 7
D1 and D2 it is observed that the characterization tests match better with the
GPR approximated design space for the lattice structures than the nonconvex
porous materials. Both Figure 7 A, B, & D indicate the importance of
additional qualitative characterizations but also prove the general accuracy
of the GPR approximation.

Although the focus of this paper is on the examination of the optimization
process for complex materials design cases, rather than simply proposing the
designs from the optimizations, we would like to still provide the eventual
extracted optimal design for each case for references. The objective values
($\mathcal{N}_{\rm bio}$), their corresponding design variables ($N_{\rm
unit}$ & $\bar{\mathcal{D}}$), and the transformed characteristic length
$\mathfrak{L}$ (in the unit of $\mu\rm m$) for all three cases benchmarked by
a nonconfined pure biofilm growth in vacuum space are shown in Table 1. Very
interestingly and unexpectedly, it is observed that all the optimal designs
extracted from porous materials confined biofilm growth exhibit more bacteria
cells in the target growth region than nonconfined biofilm growth in a vacuum
space. The optimal designs of the 2D porous membrane, lattice metamaterials,
and 3D porous media have 16%, 7%, and 11%, more biofilms in the target growth
region than the pure growth in the vacuum space, respectively. This
confinement-induced biofilm growth may help us (1) better utilize biofilms as
ELM and address the three points presented in the second paragraph in Section
1, and (2) potentially explain the natural phenomena described in the first
paragraph in Section 1. We focus on this point to conduct a further comparison
study in the following Section 3.3.

| $\mathcal{N}_{\rm bio}$ | $N_{\rm unit}$ | $\bar{\mathcal{D}}$ |  $\mathfrak{L}$ $[\mu\rm m]$
---|---|---|---|---
2D porous membrane | 32655 | 10 | 0.1 | 0.5
| 32655 | 11 | 0.1 | 0.45
Lattice metamaterials | 30096 | 1 | 0.5 | 25
3D porous media | 31152 | 7 | 1.1 | 0.71
Vacuum space | 28086 | N/A | N/A | N/A
Table 1: The highest objective values and their corresponding design variables
for different porous materials design cases, with transformed characteristic
lengths in the unit of $\mu\rm m$. For 2D porous membrane, the characteristic
length is defined as $\mathfrak{L}\equiv\mathcal{R}_{\rm vac}$. For lattice
metamaterials, the characteristic length is defined as
$\mathfrak{L}\equiv\ell_{\rm vol}$. For 3D porous media, the characteristic
length is defined as $\mathfrak{L}\equiv\mathcal{R}_{\rm vol}$. The “Vacuum
space” stand for the case where there are no porous materials defined on top
of the initial bacteria cells and their growth in a nonconfined space. Figure
7: Design space characterization for the Gaussian process regression (GPR)
reconstructed design space and topologies extraction from the characterization
processes for both the lattice metamaterials and 3D porous media design
optimization. (A1) Characterization of the design variable $N_{\rm unit}$ with
different fixed values of $\bar{\mathcal{D}}$. Note that the blue circular
dots correspond to the black triangular dots, and the grey circular dots
corresponds to the black triangular dots, in subfigure D1. The blue and red
circular dots are the characterization tests informed by qualitative
observation of the GPR reconstructed design space to approximate the optimal
design (i.e., maximal point), and the grey dots are random tests to benchmark
our characterization informed by the observations. The zoomed view describes
the detailed differences between the two sets of characterization simulations,
in which three sets of membrane topologies are selected and highlighted in red
triangular plots, nominated as $\mathbb{T}_{\alpha}$, and
$\mathbb{T}_{\beta}$, respectively. (B1) Design variable characterization for
$\bar{\mathcal{D}}$ compared with random benchmark test marked in red and grey
dots, respectively. The zoomed view describes the detailed differences between
the two sets of characterization simulations, in which three sets of membrane
topologies are selected and highlighted in red triangular plots, nominated as
$\mathbb{T}_{\beta}$, and $\mathbb{T}_{\gamma}$, respectively
($\mathbb{T}_{\beta}$ is the same topology as in subfigure A1). (C1) Extracted
porous membranes’ topologies ($\mathbb{T}_{\alpha}\sim\mathbb{T}_{\gamma}$)
from characterizing both the design variables $N_{\rm unit}$ and
$\bar{\mathcal{D}}$ corresponding to the selections in subfigures A1 & B1.
(D1) The characterization data match with the GPR reconstructed design spaces
from the EI acquisition function. The black triangular dots are the
characterization informed by observation from the GPR reconstructed design
space towards the maximal value. The grey triangular dots are randomly
selected test points to benchmark the observation-informed characterizations.
(A2) Characterization of the design variable $N_{\rm unit}$ with different
fixed values of $\bar{\mathcal{D}}$. Visualization details are the same as in
subfigure A1. (B2) Design variable characterization for $\bar{\mathcal{D}}$
compared with random benchmark test marked in red and grey dots, respectively.
Visualization details are the same as in subfigure B1, except there is no
zoomed view since the range for the objective $\mathcal{N}_{\rm bio}$ are
already within a small range. (C2) Extracted porous membranes’ topologies
($\mathbb{T}_{\alpha}\sim\mathbb{T}_{\gamma}$) from characterizing both the
design variables $N_{\rm unit}$ and $\bar{\mathcal{D}}$ corresponding to the
selections in subfigures A2 & B2. (D2) The characterization data match with
the GPR reconstructed design spaces from the EI acquisition function.
Visualization details are the same as in subfigure D1. For details please see
the texts.

###  3.3 Biomechanics of Porous Transport

Eventually, we would like to answer Question (3) in Section 1. Figure 8 shows
the benchmark study of the biofilm growth in porous membrane and vacuum space.
We pick the case of a 2D porous membrane with $N_{\rm unit}=6$ and
$\bar{\mathcal{D}}$ for comparison with biofilm growth in nonconfined vacuum
space. Figure 8 A & B visualize the snapshots of the biofilm growth
simulations, where $\tilde{\tau}$ stand for the iteration number (or time
steps), which can be converted to real-world time as $t=10\times\tilde{\tau}\
\rm[s]$. Figure 8 C visualizes the sliced view of the biofilm growth at
$\tilde{\tau}=12000$, to further explain confinement-induced biofilm growth.
Figure 8 D shows the change of the total bacteria cells $\mathcal{N}_{\rm
bio}^{\rm total}$ along with the iterations $\tilde{\tau}$, where the blue
solid line stand for biofilm growth in nonconfined vacuum space and red dashed
line stand for biofilm growth in the porous membrane. We observe two key
moments that distinguish the overall biofilm growths, the first moment is when
$\tilde{\tau}\approx 6000$ when the biofilms in the vacuum space (blue solid
line) exceeds that of in the porous materials (red dashed line); and the
second moment is when $\tilde{\tau}\approx 13500$ when the biofilms in the
porous materials (red dashed line) exceeds that of in the vacuum space (blue
solid line). The sliced views of the two moments ($\tilde{\tau}=6000$ &
$\tilde{\tau}=13500$) are visualized and indicated by shaded arrows. To
quantitatively understand the mechanism of confinement-induced biofilm growth
and transport, we compute the biofilm cell numbers distribution along the
Z-axis by counting through 100 slices at $\tilde{\tau}=12000$ (detailed
analysis can be found in ESI of Ref. Zhai and Yeo (2022)) and visualize the
results in Figure 8 E, corresponds to Figure 8 C. The blue bars indicate the
accumulative bacteria counts for biofilm growth in vacuum space and the red
bars indicate that of the porous materials.

Figure 8: Comparison study for a single 2D porous membrane with vacuum biofilm
growth case to unravel the biomechanics of porous materials induced biofilm
growth. (A) The snapshots of the simulation of biofilm growth in pure vacuum
space, where $\tilde{\tau}$ is the simulation iteration step or can be treated
as the pseudo-time. (B) The snapshots of the simulation of biofilm growth in
the 2D porous membrane. (C) Slice view of snapshot $\tilde{\tau}=12000$ for
both the 2D membrane and vacuum growth cases. (D) The accumulated bacteria
cell numbers $\mathcal{N}_{\rm bio}$ along the iteration process, where the
simulation snapshot of $\tilde{\tau}=6000$ is indicated in the left top
subfigure and $\tilde{\tau}=13500$ is indicated in the bottom right subfigure.
The solid blue line indicates the biofilm growth in vacuum space (without any
porous materials) and the red dashed line indicates the biofilm growth in the
2D porous membrane for benchmarking. The zoomed view for
$\tilde{\tau}\in[12000,15000]$ is indicated in the right subfigure with a
gradient-shaded background. (E) The bacteria cells’ spatial distribution along
the perpendicular direction (Z axis) at $\tilde{\tau}=12000$, where the cell
numbers are counted based on 100 interval slices visualized in bar plots. The
blue bars indicate the vacuum space bacteria counts and the red bars indicate
the bacteria counts in the 2D porous membrane. For details see the text.

It can be observed from Figure 8 A & B that the overall biofilms are more
densely compacted in the target growth region grew through the porous
materials compared with the vacuum space growth. From the sliced view in
Figure 8 C, we may hence propose a qualitative explanation for our
observation: the existence of the porous material takes a certain amount of
volume, which pushes the biofilm to grow upwards to occupy more space. To
break down this process more detailedly, Figure 8 D tells us that after
$\tilde{\tau}\approx 6000$ the existence of the porous materials first
suppress the biofilm growth, as $\mathcal{N}_{\rm bio}^{\rm total}$ for vacuum
space (solid blue line) first increases nonlinearly with larger values than
that of porous materials (dashed red line). But after the biofilms are well
grown in the target growth region ($\tilde{\tau}\approx 13500$), the pores in
the porous materials can be treated as “channels” that enhance the growth and
transport of biofilms. This finding is significant in the sense that the
effects of porous materials on the overall growth of biofilms change in
different stages of the growth processes within the pores. Based on these
comprehensive qualitative analyses, Figure 8 E offers quantitative evidence
that porous material push biofilms’ upward growth by taking up volumetric
spaces — the biofilm accumulation within the porous materials spatial range
($Z\in[12.5,37.5]\mu\rm m$) for porous materials (red bars) are evidently
smaller than that of vacuum space (blue bars). Based on the bacteria cell
numbers count from 100 slices, the porous region bacteria counts for porous
membrane and vacuum space are 48643 and 58482, respectively, where the vacuum
space contains 20% more biofilms than that constrained in the porous membrane.
The target growth region bacteria counts for porous membrane and vacuum space
are 31404 and 13764, respectively, where the porous membrane contains 128%
more biofilms than that grew in the vacuum space. The data not only verifies
our qualitative explanations that the porous membrane facilitates biofilm
growth by taking up volumetric space but also further explains how the porous
membrane increases the overall biofilms — the pores behave like channels that
transport biofilms to the target region so that the bacteria count in the
target growth region for porous membrane are significantly large than that of
in vacuum space.

##  4 Conclusions & Outlook

In this paper, we present efforts to design different porous materials for
enhanced biofilm transport and control from computational models using
Bayesian optimization. We focus on characterizing the design optimization
process, comprehensively analyzing the approximated design space, and further
providing in-depth physical insights from the optimization. We formulate three
different types of porous structural materials for design optimization aiming
to maximize the biofilms in the target growth region. For three different
types of porous materials, the trends of the reconstructed design space match
well with the sampling density. For the 2D porous membrane, the variance of
the overall samples by the UCB acquisition function is 32.08% higher than that
of the EI acquisition function; the mean objective of the overall samples by
the EI acquisition function is 1.49% higher than that of the UCB acquisition
function. Given the predefined target region of higher sampled densities, the
EI acquisition function is 2.35% more efficient than the UCB acquisition
function compared with uniformly distributed grid search methods by estimating
the last 100 sampling points. The GPR approximated design spaces match well
with the selected characterization tests. Using only the EI acquisition
function, we conduct the design space characterization for lattice
metamaterials and porous media under the same procedure. For the lattice
metamaterials, by looking at the last 100 samples in the predefined target
design space, BO is 92.89% more efficient than the uniform grid search. For
the 3D porous media, there are 223.04% more sampled points by BO than the
uniform grid search in the predefined target design space. We further provide
the design variables of the selected optimal design for different porous
materials formulations. Very interestingly, all the extracted optimal designs
have more bacteria cells in the target growth region than pure biofilm growth
in the vacuum space without any confinement. We conduct a comparison study
trying to understand the phenomenon and found that there are 20% more biofilms
for the vacuum space than that confined in the porous materials. What’s more,
there are 128% more biofilms in the target growth region for the porous
materials-induced biofilm growth compared with the vacuum space growth. We
thence propose that the existence of porous materials stimulates the biofilms
by taking up volumetric space to push upward growth. Note that this is not
universally tested for all kinds of porous materials with all radii range, and
testing the size effect for confinement-induced biofilm growth would be our
follow-up work in the future.

Our work is significant and innovative from three major aspects: (1)
Implications and guidance to broad audiences. Our work could inspire theorists
and programmers to develop new theories and algorithms for modeling biofilm
and guide experimentalists to conduct new investigations. (2) Rigorous and
comprehensive optimization analysis of the optimization process and direct
characterization of the design space. (3) Understanding the mechanism from
both the optimization characterization and computational modeling brings in
new knowledge. From both three aspects, our work reaches a broad range of
different research areas spanning mechanics, materials, machine learning,
biology, environments, and many fields. This paper, to our knowledge, is the
first work that utilizes ML as an optimization tool for characterizing the
underlying mechanisms of confined biofilm dynamics using computational models.
Our work is expected to unveil a new paradigm of conducting inverse design to
inspire physics discovery by leveraging computational models, ML, and design
optimizations.

## Acknowledgement

J.Y. acknowledges support from the US National Science Foundation (grant nos.
2038057 and 2223785) and the Cornell University faculty startup grant. The
authors also acknowledge the computational resources provided by the XSEDE
program under grants TG-MAT200004 and TG-BIO210063 and the computational
resources provided by the G2 cluster from Cornell University.

## Appendix

Supplementary Figure 1 visualizes the overall design processes for lattice
metamaterials and 3D porous media, respectively (Figure 2 B & C). The upper
figure (1) stands for the change of the objectives and the lower figure stand
for the design variables’ changes w.r.t. the iterations, similar to what has
been shown in Figure 3. A converging process of the objective values is
observed for 3D porous media (B1), whereas the objectives are most fluctuating
more for the lattice metamaterials (A1), which can be attributed to the
nonconvex design space in Figure 6. For the lattice metamaterials, the design
variables are fluctuating along the iterations where $\bar{\mathcal{D}}$ is
sampled toward higher values and $N_{\rm unit}$ is sampled toward the lower
(A2). For the 3D porous media, similar trends are also observed yet the
difference is they are initially sampled in a similar value range and the
discrepancy of the sampling value trends begins to occur after approximately
300 iterations (B2).

Supplementary Figure 1: Design optimization process for both the lattice
metamaterials and 3D porous media. The top subfigures indicate the objective
value (bacteria cell numbers $\mathcal{N}_{\rm bio}$) along the simulation
iteration processes. The bottom subfigures indicate the normalized design
variable maps to visualize the changing trends of the design variables. (A1)
Objective value change for the lattice metamaterials optimization case. (A2)
Objective value change for the 3D porous media optimization case. (B1) The
design variables’ change for the lattice metamaterials optimization case. (B2)
The design variables’ change for the 3D porous media optimization case.

## References

  * Costerton et al. [1999] J. W. Costerton, Philip S. Stewart, and E. P. Greenberg.  Bacterial biofilms: A common cause of persistent infections.  _Science_ , 284(5418):1318–1322, May 1999.  doi:10.1126/science.284.5418.1318.  URL https://doi.org/10.1126/science.284.5418.1318.
  * Friedlander et al. [2013] Ronn S. Friedlander, Hera Vlamakis, Philseok Kim, Mughees Khan, Roberto Kolter, and Joanna Aizenberg.  Bacterial flagella explore microscale hummocks and hollows to increase adhesion.  _Proceedings of the National Academy of Sciences_ , 110(14):5624–5629, March 2013.  doi:10.1073/pnas.1219662110.  URL https://doi.org/10.1073/pnas.1219662110.
  * Yebra et al. [2004] Diego Meseguer Yebra, Søren Kiil, and Kim Dam-Johansen.  Antifouling technology—past, present and future steps towards efficient and environmentally friendly antifouling coatings.  _Progress in Organic Coatings_ , 50(2):75–104, July 2004.  doi:10.1016/j.porgcoat.2003.06.001.  URL https://doi.org/10.1016/j.porgcoat.2003.06.001.
  * Dobretsov et al. [2006] Sergey Dobretsov, Hans-Uwe Dahms, and PerI-Yuan Qian.  Inhibition of biofouling by marine microorganisms and their metabolites.  _Biofouling_ , 22(1):43–54, January 2006.  doi:10.1080/08927010500504784.  URL https://doi.org/10.1080/08927010500504784.
  * Donlan and Costerton [2002] Rodney M. Donlan and J. William Costerton.  Biofilms: Survival mechanisms of clinically relevant microorganisms.  _Clinical Microbiology Reviews_ , 15(2):167–193, April 2002.  doi:10.1128/cmr.15.2.167-193.2002.  URL https://doi.org/10.1128/cmr.15.2.167-193.2002.
  * Jonkers et al. [2010] Henk M. Jonkers, Arjan Thijssen, Gerard Muyzer, Oguzhan Copuroglu, and Erik Schlangen.  Application of bacteria as self-healing agent for the development of sustainable concrete.  _Ecological Engineering_ , 36(2):230–235, February 2010.  doi:10.1016/j.ecoleng.2008.12.036.  URL https://doi.org/10.1016/j.ecoleng.2008.12.036.
  * Chattopadhyay et al. [2022] Indranil Chattopadhyay, Rajesh Banu J, T. M. Mohamed Usman, and Sunita Varjani.  Exploring the role of microbial biofilm for industrial effluents treatment.  _Bioengineered_ , 13(3):6420–6440, February 2022\.  doi:10.1080/21655979.2022.2044250.  URL https://doi.org/10.1080/21655979.2022.2044250.
  * Balasubramanian et al. [2019] Srikkanth Balasubramanian, Marie-Eve Aubin-Tam, and Anne S. Meyer.  3d printing for the fabrication of biofilm-based functional living materials.  _ACS Synthetic Biology_ , 8(7):1564–1567, July 2019.  doi:10.1021/acssynbio.9b00192.  URL https://doi.org/10.1021/acssynbio.9b00192.
  * Zhai and Yeo [2022] Hanfeng Zhai and Jingjie Yeo.  Computational design of antimicrobial active surfaces via automated bayesian optimization.  _ACS Biomaterials Science & Engineering_, 9(1):269–279, Dec 2022.  ISSN 2373-9878.  doi:10.1021/acsbiomaterials.2c01079.  URL http://dx.doi.org/10.1021/acsbiomaterials.2c01079.
  * Liu et al. [2022] Xiaomeng Liu, Toshiyuki Ueki, Hongyan Gao, Trevor L. Woodard, Kelly P. Nevin, Tianda Fu, Shuai Fu, Lu Sun, Derek R. Lovley, and Jun Yao.  Microbial biofilms for electricity generation from water evaporation and power to wearables.  _Nature Communications_ , 13(1), July 2022.  doi:10.1038/s41467-022-32105-6.  URL https://doi.org/10.1038/s41467-022-32105-6.
  * Mukherjee and Cao [2020] Manisha Mukherjee and Bin Cao.  Engineering controllable biofilms for biotechnological applications.  _Microbial Biotechnology_ , 14(1):74–78, November 2020.  doi:10.1111/1751-7915.13715.  URL https://doi.org/10.1111/1751-7915.13715.
  * Powell et al. [2018] Lydia C. Powell, Manon F. Pritchard, Elaine L. Ferguson, Kate A. Powell, Shree U. Patel, Phil D. Rye, Stavroula-Melina Sakellakou, Niklaas J. Buurma, Charles D. Brilliant, Jack M. Copping, Georgina E. Menzies, Paul D. Lewis, Katja E. Hill, and David W. Thomas.  Targeted disruption of the extracellular polymeric network of pseudomonas aeruginosa biofilms by alginate oligosaccharides.  _npj Biofilms and Microbiomes_ , 4(1), June 2018.  doi:10.1038/s41522-018-0056-3.  URL https://doi.org/10.1038/s41522-018-0056-3.
  * Xu et al. [2011] Zhijie Xu, Paul Meakin, Alexandre Tartakovsky, and Timothy D. Scheibe.  Dissipative-particle-dynamics model of biofilm growth.  _Physical Review E_ , 83(6), June 2011.  doi:10.1103/physreve.83.066702.  URL https://doi.org/10.1103/physreve.83.066702.
  * Brandani et al. [2015] Giovanni B. Brandani, Marieke Schor, Ryan Morris, Nicola Stanley-Wall, Cait E. MacPhee, Davide Marenduzzo, and Ulrich Zachariae.  The bacterial hydrophobin BslA is a switchable ellipsoidal janus nanocolloid.  _Langmuir_ , 31(42):11558–11563, October 2015\.  doi:10.1021/acs.langmuir.5b02347.  URL https://doi.org/10.1021/acs.langmuir.5b02347.
  * Smith et al. [2007] Bryan Smith, Benjamin Vaughan, and David Chopp.  The extended finite element method for boundary layer problems in biofilm growth.  _Communications in Applied Mathematics and Computational Science_ , 2(1):35–56, August 2007.  doi:10.2140/camcos.2007.2.35.  URL https://doi.org/10.2140/camcos.2007.2.35.
  * Li et al. [2019] Bowen Li, Denis Taniguchi, Jayathilake Pahala Gedara, Valentina Gogulancea, Rebeca Gonzalez-Cabaleiro, Jinju Chen, Andrew Stephen McGough, Irina Dana Ofiteru, Thomas P. Curtis, and Paolo Zuliani.  NUFEB: A massively parallel simulator for individual-based modelling of microbial communities.  _PLOS Computational Biology_ , 15(12):e1007125, December 2019.  doi:10.1371/journal.pcbi.1007125.  URL https://doi.org/10.1371/journal.pcbi.1007125.
  * Kapellos et al. [2015] George E. Kapellos, Terpsichori S. Alexiou, and Stavros Pavlou.  Fluid-biofilm interactions in porous media.  In _Heat Transfer and Fluid Flow in Biological Processes_ , pages 207–238. Elsevier, 2015.  doi:10.1016/b978-0-12-408077-5.00008-0.  URL https://doi.org/10.1016/b978-0-12-408077-5.00008-0.
  * Galy et al. [2012] Olivier Galy, Patricia Latour-Lambert, Kais Zrelli, Jean-Marc Ghigo, Christophe Beloin, and Nelly Henry.  Mapping of bacterial biofilm local mechanics by magnetic microparticle actuation.  _Biophysical Journal_ , 103(6):1400–1408, September 2012.  doi:10.1016/j.bpj.2012.07.001.  URL https://doi.org/10.1016/j.bpj.2012.07.001.
  * Rodrigo-Navarro et al. [2021] Aleixandre Rodrigo-Navarro, Shrikrishnan Sankaran, Matthew J. Dalby, Aránzazu del Campo, and Manuel Salmeron-Sanchez.  Engineered living biomaterials.  _Nature Reviews Materials_ , 6(12):1175–1190, August 2021.  doi:10.1038/s41578-021-00350-8.  URL https://doi.org/10.1038/s41578-021-00350-8.
  * Li et al. [2021] Nan K. Li, Yuxin Xie, and Yaroslava G. Yingling.  Insights into structure and aggregation behavior of elastin-like polypeptide coacervates: All-atom molecular dynamics simulations.  _The Journal of Physical Chemistry B_ , 125(30):8627–8635, July 2021.  doi:10.1021/acs.jpcb.1c02822.  URL https://doi.org/10.1021/acs.jpcb.1c02822.
  * Duddu et al. [2008] Ravindra Duddu, Stéphane Bordas, David Chopp, and Brian Moran.  A combined extended finite element and level set method for biofilm growth.  _International Journal for Numerical Methods in Engineering_ , 74(5):848–870, 2008.  doi:10.1002/nme.2200.  URL https://doi.org/10.1002/nme.2200.
  * Zhai et al. [2022] Hanfeng Zhai, Quan Zhou, and Guohui Hu.  Predicting micro-bubble dynamics with semi-physics-informed deep learning.  _AIP Advances_ , 12(3), 03 2022.  ISSN 2158-3226.  doi:10.1063/5.0079602.  URL https://doi.org/10.1063/5.0079602.  035153\.
  * Hadamard [1902] J. Hadamard.  On problems with partial derivatives and their physical significance.  _Princeton University Bulletin_ , 13:49–52, 1902.
  * Frazier [2018] Peter I. Frazier.  A tutorial on bayesian optimization, 2018.  URL https://arxiv.org/abs/1807.02811.
  * Fuhg and Bouklas [2022] Jan N. Fuhg and Nikolaos Bouklas.  On physics-informed data-driven isotropic and anisotropic constitutive models through probabilistic machine learning and space-filling sampling.  _Computer Methods in Applied Mechanics and Engineering_ , 394:114915, May 2022.  doi:10.1016/j.cma.2022.114915.  URL https://doi.org/10.1016/j.cma.2022.114915.
  * Sutton and Barto [2018] Richard S Sutton and Andrew G Barto.  _Reinforcement Learning: An Introduction_.  MIT Press, 2018.  URL https://web.stanford.edu/class/psych209/Readings/SuttonBartoIPRLBook2ndEd.pdf.
  * Mitchell [1998] Melanie Mitchell.  _An Introduction to Genetic Algorithms_.  MIT Press, 1998.
  * Kennedy and Eberhart [1995] J. Kennedy and R. Eberhart.  Particle swarm optimization.  In _Proceedings of ICNN’95 - International Conference on Neural Networks_ , volume 4, pages 1942–1948 vol.4, 1995.  doi:10.1109/ICNN.1995.488968.
  * Monod [1949] Jacques Monod.  THE GROWTH OF BACTERIAL CULTURES.  _Annual Review of Microbiology_ , 3(1):371–394, October 1949.  doi:10.1146/annurev.mi.03.100149.002103.  URL https://doi.org/10.1146/annurev.mi.03.100149.002103.
  * Xavier et al. [2005] Joao B. Xavier, Cristian Picioreanu, and Mark C. M. van Loosdrecht.  A framework for multidimensional modelling of activity and structure of multispecies biofilms.  _Environmental Microbiology_ , 7(8):1085–1103, August 2005.  doi:10.1111/j.1462-2920.2005.00787.x.  URL https://doi.org/10.1111/j.1462-2920.2005.00787.x.
  * Jayathilake et al. [2017] Pahala Gedara Jayathilake, Prashant Gupta, Bowen Li, Curtis Madsen, Oluwole Oyebamiji, Rebeca González-Cabaleiro, Steve Rushton, Ben Bridgens, David Swailes, Ben Allen, A. Stephen McGough, Paolo Zuliani, Irina Dana Ofiteru, Darren Wilkinson, Jinju Chen, and Tom Curtis.  A mechanistic individual-based model of microbial communities.  _PLOS ONE_ , 12(8):e0181965, August 2017.  doi:10.1371/journal.pone.0181965.  URL https://doi.org/10.1371/journal.pone.0181965.
  * Biswas et al. [2021] Arpan Biswas, Anna N. Morozovska, Maxim Ziatdinov, Eugene A. Eliseev, and Sergei V. Kalinin.  Multi-objective bayesian optimization of ferroelectric materials with interfacial control for memory and energy storage applications.  _Journal of Applied Physics_ , 130(20):204102, November 2021.  doi:10.1063/5.0068903.  URL https://doi.org/10.1063/5.0068903.
  * Deshwal et al. [2021] Aryan Deshwal, Cory M. Simon, and Janardhan Rao Doppa.  Bayesian optimization of nanoporous materials.  _Molecular Systems Design & Engineering_, 6(12):1066–1086, 2021.  doi:10.1039/d1me00093d.  URL https://doi.org/10.1039/d1me00093d.
  * Snoek et al. [2012] Jasper Snoek, Hugo Larochelle, and Ryan P. Adams.  Practical bayesian optimization of machine learning algorithms, 2012.  URL https://arxiv.org/abs/1206.2944.
  * Pousti et al. [2019] Mohammad Pousti, Mir Pouyan Zarabadi, Mehran Abbaszadeh Amirdehi, François Paquet-Mercier, and Jesse Greener.  Microfluidic bioanalytical flow cells for biofilm studies: a review.  _The Analyst_ , 144(1):68–86, 2019.  doi:10.1039/c8an01526k.  URL https://doi.org/10.1039/c8an01526k.
  * Ye et al. [2021] Dingding Ye, Pengqing Zhang, Jun Li, Xun Zhu, Rong Chen, and Qiang Liao.  In situ visualization of biofilm formation in a microchannel for a microfluidic microbial fuel cell anode.  _International Journal of Hydrogen Energy_ , 46(27):14651–14658, April 2021.  doi:10.1016/j.ijhydene.2020.08.170.  URL https://doi.org/10.1016/j.ijhydene.2020.08.170.
  * Landa-Marbán et al. [2019] David Landa-Marbán, Iuliu Sorin Pop, Kundan Kumar, and Florin A. Radu.  Numerical simulation of biofilm formation in a microchannel.  In _Lecture Notes in Computational Science and Engineering_ , pages 799–807. Springer International Publishing, 2019.  doi:10.1007/978-3-319-96415-7_75.  URL https://doi.org/10.1007/978-3-319-96415-7_75.
  * Aspa et al. [2011] Y. Aspa, G. Debenest, and M. Quintard.  Effective dispersion in channelled biofilms.  _International Journal of Environment and Waste Management_ , 7(1/2):112, 2011.  doi:10.1504/ijewm.2011.037371.  URL https://doi.org/10.1504/ijewm.2011.037371.
  * Landa-Marbán et al. [2020] David Landa-Marbán, Gunhild Bødtker, Kundan Kumar, Iuliu S. Pop, and Florin A. Radu.  An upscaled model for permeable biofilm in a thin channel and tube.  _Transport in Porous Media_ , 132(1):83–112, January 2020.  doi:10.1007/s11242-020-01381-5.  URL https://doi.org/10.1007/s11242-020-01381-5.
  * Ma et al. [2022] Chunping Ma, Yilong Chang, Shuai Wu, and Ruike Renee Zhao.  Deep learning-accelerated designs of tunable magneto-mechanical metamaterials.  _ACS Applied Materials & Interfaces_, 14(29):33892–33902, July 2022.  doi:10.1021/acsami.2c09052.  URL https://doi.org/10.1021/acsami.2c09052.
  * Shaw et al. [2019] Lucas A. Shaw, Frederick Sun, Carlos M. Portela, Rodolfo I. Barranco, Julia R. Greer, and Jonathan B. Hopkins.  Computationally efficient design of directionally compliant metamaterials.  _Nature Communications_ , 10(1), January 2019.  doi:10.1038/s41467-018-08049-1.  URL https://doi.org/10.1038/s41467-018-08049-1.
  * Gu [2018] X. Wendy Gu.  Mechanical properties of architected nanomaterials made from organic–inorganic nanocrystals.  _JOM_ , 70(10):2205–2217, August 2018.  doi:10.1007/s11837-018-3094-7.  URL https://doi.org/10.1007/s11837-018-3094-7.
  * Portela et al. [2020] Carlos M. Portela, A. Vidyasagar, Sebastian Krödel, Tamara Weissenbach, Daryl W. Yee, Julia R. Greer, and Dennis M. Kochmann.  Extreme mechanical resilience of self-assembled nanolabyrinthine materials.  _Proceedings of the National Academy of Sciences_ , 117(11):5686–5693, March 2020.  doi:10.1073/pnas.1916817117.  URL https://doi.org/10.1073/pnas.1916817117.
  * Ovelheiro [2020] Bryan Ovelheiro.  3d printed architected materials for improving biofilm carriers for wastewater treatment applications.  Masters project, University of Massachusetts Amherst, Amherst, MA, 2020\.  URL https://scholarworks.umass.edu/cee_masters/20/.
  * He et al. [2021] Yu-Ting He, Qian Fu, Yuan Pang, Qing Li, Jun Li, Xun Zhu, Ren-Hao Lu, Wei Sun, Qiang Liao, and Uwe Schröder.  Customizable design strategies for high-performance bioanodes in bioelectrochemical systems.  _iScience_ , 24(3):102163, March 2021.  doi:10.1016/j.isci.2021.102163.  URL https://doi.org/10.1016/j.isci.2021.102163.
  * Hall et al. [2021] Donald C. Hall, Phillip Palmer, Hai-Feng Ji, Garth D. Ehrlich, and Jarosław E. Król.  Bacterial biofilm growth on 3d-printed materials.  _Frontiers in Microbiology_ , 12, May 2021.  doi:10.3389/fmicb.2021.646303.  URL https://doi.org/10.3389/fmicb.2021.646303.
  * Bhattacharjee and Datta [2019] Tapomoy Bhattacharjee and Sujit S. Datta.  Bacterial hopping and trapping in porous media.  _Nature Communications_ , 10(1), May 2019.  doi:10.1038/s41467-019-10115-1.  URL https://doi.org/10.1038/s41467-019-10115-1.
  * Carrel et al. [2018] Maxence Carrel, Verónica L. Morales, Mario A. Beltran, Nicolas Derlon, Rolf Kaufmann, Eberhard Morgenroth, and Markus Holzner.  Biofilms in 3d porous media: Delineating the influence of the pore network geometry, flow and mass transfer on biofilm development.  _Water Research_ , 134:280–291, May 2018.  doi:10.1016/j.watres.2018.01.059.  URL https://doi.org/10.1016/j.watres.2018.01.059.
  * Coyte et al. [2016] Katharine Z. Coyte, Hervé Tabuteau, Eamonn A. Gaffney, Kevin R. Foster, and William M. Durham.  Microbial competition in porous environments can select against rapid biofilm growth.  _Proceedings of the National Academy of Sciences_ , 114(2), December 2016.  doi:10.1073/pnas.1525228113.  URL https://doi.org/10.1073/pnas.1525228113.
  * Kurz et al. [2022] Dorothee L. Kurz, Eleonora Secchi, Francisco J. Carrillo, Ian C. Bourg, Roman Stocker, and Joaquin Jimenez-Martinez.  Competition between growth and shear stress drives intermittency in preferential flow paths in porous medium biofilms.  _Proceedings of the National Academy of Sciences_ , 119(30), July 2022.  doi:10.1073/pnas.2122202119.  URL https://doi.org/10.1073/pnas.2122202119.
  * Dehkharghani et al. [2023] Amin Dehkharghani, Nicolas Waisbord, and Jeffrey S. Guasto.  Self-transport of swimming bacteria is impaired by porous microstructure.  _Communications Physics_ , 6(1), January 2023.  doi:10.1038/s42005-023-01136-w.  URL https://doi.org/10.1038/s42005-023-01136-w.