# VoroTO: Multiscale Topology Optimization of Voronoi Structures using
Surrogate Neural Networks

Rahul Kumar Padhy
Department of Mechanical Engineering
University of Wisconsin-Madison
Madison, WI, USA
<EMAIL_ADDRESS>
Krishnan Suresh
Department of Mechanical Engineering
University of Wisconsin-Madison
Madison, WI, USA
<EMAIL_ADDRESS>
Aaditya Chandrasekhar
Department of Mechanical Engineering
University of Wisconsin-Madison
Madison, WI, USA
<EMAIL_ADDRESS>

###### Abstract

Cellular structures found in nature exhibit remarkable properties such as high
strength, high energy absorption, excellent thermal/acoustic insulation, and
fluid transfusion. Many of these structures are Voronoi-like; therefore
researchers have proposed Voronoi multi-scale designs for a wide variety of
engineering applications. However, designing such structures can be
computationally prohibitive due to the multi-scale nature of the underlying
analysis and optimization. In this work, we propose the use of a neural
network (NN) to carry out efficient topology optimization (TO) of multi-scale
Voronoi structures. The NN is first trained using Voronoi parameters (cell
site locations, thickness, orientation, and anisotropy) to predict the
homogenized constitutive properties. This network is then integrated into a
conventional TO framework to minimize structural compliance subject to a
volume constraint. Special considerations are given for ensuring positive
definiteness of the constitutive matrix and promoting macroscale connectivity.
Several numerical examples are provided to showcase the proposed method.

_K_ eywords Topology Optimization $\cdot$ Voronoi Structures $\cdot$ Neural
Network

Figure 1: Graphical abstract: Offline computation: Given a dataset containing
Voronoi microstructure parameters and homogenized constitutive properties, a
neural network is trained offline. Multiscale TO: The trained network is used
as a surrogate during topology optimization to derive optimized Voronoi
structures.

##  1 Introduction

Topology optimization is a powerful engineering tool employed to optimize the
distribution of material within a given design domain to achieve optimal
performance under specified constraints [1, 2]. Traditionally, single-scale
topology optimization focuses on determining the optimal layout of material
within a structure at a single scale, aiming primarily at maximizing stiffness
[1] for a given set of loading conditions. However, engineering applications
often demand more than just structural efficiency [3]. Features such as
structural robustness [4], high strength [5], superior energy absorption [6],
fluid circulation [7], thermal and acoustic insulation [8] capabilities are
crucial. These attributes not only enhance the overall performance of
structures but also facilitate multifunctionality, making them indispensable
across various industries, ranging from aerospace [9, 10], biomedical [11,
12], and more [13].

Many natural structures, such as bones [14], wood [15], and insect wings [16]
exhibit remarkable combinations of strength, energy absorption, fluid
circulation, and insulating properties. These features are often achieved
through the intricate arrangement of porous structures, which possess
desirable properties such as stochasticity [17], anisotropy [18, 19], and
connectivity [20]. For instance, consider the calcaneus (heel bone) and its
associated loading conditions during walking, as depicted in fig. 2. This
scenario can be effectively modeled as a simple cantilever beam as in fig. 3a
[19]. The growth of bone entails a continuous self-optimization process, where
its internal structure adapts to maximize mechanical efficiency under varying
external conditions. This behavior can be modeled through topology
optimization, by formulating the problem of maximizing stiffness [21, 19].

Figure 2: Heel bone and loading conditions.

The heel bone topology optimization problem is idealized fig. 3a. An optimal
design that maximizes stiffness using single-scale topology optimization is
presented in fig. 3b. Although this design is optimal for stiffness, it lacks
the porosity essential for fluid movement. In bones, pores play a vital role
in blood and interstitial fluid movement, enabling the exchange of oxygen and
nutrients between tissues [7]. A hypothetical porous structure that sacrifices
stiffness for porosity [22, 23, 24] is depicted in fig. 3c.

Systematic computation of such structures poses significant computational
challenges due to their multiscale nature [25, 26]. Theoretically, one can use
high-resolution single-scale topology optimization to arrive at such
structures, but this is not viable in practice [26]. Here, we present a neural
network (NN) based computationally efficient multiscale topology optimization
approach for computing such structures. Furthermore, we use Voronoi
microstructures to represent porosity as they offer high design freedom [27],
exhibit anisotropy, and often resemble the porous structures found in bone
[28, 29]. A parametric representation of the microstructures is used to train
an NN, facilitating the mapping of microstructure parameters to homogenized
properties. Special attention is given to the parametric representation to
promote macroscale connectivity, addressing a limitation [30] not inherently
present in unit cell-based multiscale topology optimization approaches [31,
32].

Figure 3: (a) Topology optimization problem. (b) Single scale design. (c) A
multiscale porous design.

###  1.1 Topology optimization of porous structures

Topology optimization of porous structures has gained considerable attention,
with two primary approaches [30]: full-scale and multi-scale. The full-scale
approach employs geometrical constraints such as pattern repetition [33, 34]
and local volume constraints [35, 36, 37] to generate optimized porous
structures. Recent enhancements include incorporating length scale constraint
[38] and parametrized porous structures [39]. Although these approaches can
generate optimized porous structures, the incorporation of additional
geometrical constraints often results in reduced performance when compared to
traditional topology optimization methods [40]. Furthermore, these approaches
can be computationally expensive because of the need for high-resolution
topology optimization of fine porous structures [26].

Multi-scale approaches separate the optimization problem into macro and micro
scales, enabling customization of local material properties and macroscale
structural behavior [30, 41]. It has yielded materials with exceptional
properties, such as micro-modules with high Young’s modulus [42]. Advanced
techniques such as clustering [43] and level-set-based [44] optimization are
used to enhance structural properties and to obtain designs with smooth
boundaries respectively. Fixed types of parametric cells were used to design
porous structures such as planar rod networks [45], or TPMS (triply periodic
minimal surfaces) based [46, 47, 48]. These approaches reduce computational
costs in generating foams with valid geometry but have a limited design space
[49]. For a comprehensive review of these approaches, please see [30].

###  1.2 Voronoi based porous structure optimization

Traditionally, regular shapes such as diamond [50] and honeycomb [51] have
been used to represent porous structures. However, these regular structures
are susceptible to stress shielding and often exhibit limited permeability for
fluid movement [52]. Consequently, there is growing research interest in
irregular porous structures, specifically, Voronoi tessellation [53, 27], to
model the complex shapes of bone structures. Voronoi tessellations offer high
design freedom [27], exhibit anisotropy, and often resemble the porous
structures found in nature [28, 29].

Various stochastic methods based on Voronoi tessellations have been proposed
[54, 55, 26], but these lead to extensive data occupancy and are
computationally intensive [56]. An alternative method involves utilizing seed
positions and beam radii as design parameters, enabling them to conform to
freeform shapes while maintaining geometric connectivity. However, this does
not facilitate the generation of anisotropic designs [49]. Another technique
employs an iterative process of removing and adding site points to achieve a
desired stress distribution on the porous structure, but it lacks the
capability for optimizing a Voronoi structure through continuous sensitivity
analysis [57]. A multi-scale approach was introduced for designing Voronoi-
graded cellular structures for heat transfer problems. This method adapts
offline computation to reduce computational expenses [58]. Nevertheless, the
resulting design lacks anisotropy and has a constant Voronoi strut thickness
across the entire design domain.

###  1.3 Contributions

This paper presents a novel framework for generating optimized Voronoi
multiscale designs with spatially varying thickness and anisotropy. In
contrast to full-scale approach for Voronoi design generation [39], we
leverage a multiscale design paradigm, utilizing an offline-trained NN that
maps Voronoi microstructure parameters to homogenized constitutive
microstructure properties. Key contributions of this work include:

  * •

Design space: We consider anisotropic Voronoi cells of varying thickness in
our design space, enhancing performance [29].

  * •

Connectivity: In the process of designing Voronoi microstructures, neighboring
cell sites are taken into account to promote connectivity.

  * •

Physically valid: We ensure that the stiffness matrix predicted by the network
is positive definite, i.e., physically valid.

  * •

End-to-end differentiability: Our methodology is end-to-end differentiable,
enabling automated sensitivity analysis for gradient-based optimization.

##  2 Proposed Method

###  2.1 Voronoi Diagram and Structure

Our objective here is to design porous structures using Voronoi diagrams. A
Voronoi diagram is a geometric structure that partitions space into regions or
cells based on the proximity to a specified set of points called sites [59,
60, 61]; see fig. 4. The shape of the Voronoi cell associated with a site is
only influenced by that site and nearby sites. Voronoi diagrams are valuable
for representing and analyzing complex spatial relationships, making them
essential in computational geometry [62, 63], modeling of porous structures
[53], and design problems [39].

Figure 4:  Voronoi diagram defined by a set of sites (points). The shape of
cell ($\bigstar$) is influenced by its immediate neighbors
($\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}$).

Using the Voronoi _diagram_ as a basis, we will now define a Voronoi
_structure_. Let the domain contain $\tilde{S}$ sites (points), denoted by
$(\bar{x}_{s},\bar{y}_{s});s=1,2,...,\tilde{S}$; this results in $\tilde{S}$
cells. We now associate 3 parameters with each cell: (1) $\beta_{s}$ to create
a thickness to the cells, (2) $\alpha_{s}$ to control the degree of
anisotropy, and (3) $\theta_{s}$ to control the orientation. These parameters
are essential to obtain an optimal distribution of porous material. The
resulting Voronoi structure is defined using a density function as follows.
Given any point $(x,y)$ we define an anisotropic [64, 65] and oriented
distance [66] of $(x,y)$ to any cell site $s$ as:

$d_{s}(x,y)=\sqrt{\alpha_{s}(\delta_{s}^{x})^{2}+(\delta_{s}^{y})^{2}/\alpha_{s}}$
(1)

where,

$\begin{pmatrix}\delta_{s}^{x}\\\
\delta_{s}^{y}\end{pmatrix}=\begin{bmatrix}\cos\theta_{s}&-\sin\theta_{s}\\\
\sin\theta_{s}&\cos\theta_{s}\end{bmatrix}\begin{pmatrix}x-\bar{x}_{s}\\\
y-\bar{y}_{s}\end{pmatrix}$ (2)

Then the density field associated with a site $s$ is defined via a soft-max
function [39]:

$\hat{\rho}_{s}(x,y)=\bigg{(}\frac{e^{-kd_{s}(x,y)}}{\sum\limits_{s=1}^{S}e^{-kd_{s}(x,y)}}\bigg{)}^{\beta_{s}}$
(3)

where the parameter $k$ controls the sharpness of the soft-max function [39].
Observe that $\hat{\rho}_{s}(\bar{x}_{s},\bar{y}_{s})=1$ when the point ($x$,
$y$) is closest to site $s$, and approximately 0 far away from it. Finally, we
can compute the total density field at ($x$, $y$) as:

$\rho(x,y)=1-\sum\limits_{s=1}^{\tilde{S}}\hat{\rho}_{s}(x,y)$ (4)

A typical density field is illustrated in fig. 5. While the sites control the
topology of the Voronoi structure, the parameters $\beta_{s}$, $\alpha_{s}$,
and $\theta_{s}$ control its shape; additional examples are provided later.

Figure 5: A typical density field defining a Voronoi structure.

###  2.2 Single Scale Voronoi Structure Optimization

In principle, one can now optimize the location of the $\tilde{S}$ cell sites
$(\bar{x}_{s},\bar{y}_{s})$, the thicknesses $\beta_{s}$, anisotropies
$\alpha_{s}$ and orientations $\theta_{s}$ to design optimal Voronoi
structures. For example, one may pose a compliance minimization problem,
subject to a volume constraint as:

$\displaystyle\underset{\bar{\bm{x}},\bar{\bm{y}},\bm{\beta},\bm{\theta},\bm{\alpha}}{\text{minimize}}$
$\displaystyle J=\bm{f}^{T}\bm{u}$ (5a) subject to
$\displaystyle\bm{K}\bm{u}=\bm{f}$ (5b) and $\displaystyle V\leq V^{max}$ (5c)
$\displaystyle x_{min}\leq\bar{x}_{s}\leq x_{max}\;,\;\forall s$ (5d)
$\displaystyle y_{min}\leq\bar{y}_{s}\leq y_{max}\;,\;\forall s$ (5e)
$\displaystyle\beta_{min}\leq\beta_{s}\leq\beta_{max}\;,\;\forall s$ (5f)
$\displaystyle\alpha_{min}\leq\alpha_{s}\leq\alpha_{max}\;,\;\forall s$ (5g)
$\displaystyle\theta_{min}\leq\theta_{s}\leq\theta_{max}\;,\;\forall s$ (5h)

where $\bm{K}$ is global stiffness matrix, $\bm{f}$ is the applied load. For
example, the authors of [39] illustrate optimizing Voronoi structures by
varying cell site locations, thickness, and anisotropy to minimize compliance.
However, this is computationally prohibitive for designing fine-scale porous
structures since it will require a large number of cell sites $\tilde{S}$,
but, more importantly, it will entail a dense finite element mesh to capture
the thin features. In other words, a single-scale optimization is not a viable
strategy for achieving porous structures. Instead, we will pursue a multi-
scale framework, discussed in the next section.

###  2.3 Proposed Method: Overview

To achieve a multi-scale design, we discretize the domain into a finite number
of _macro_ elements, where each element $e$ contains a small number of $S$
cell sites, $(\bar{x}^{(e)}_{s},\bar{y}^{(e)}_{s})$; $s=1,2,...,S$; $S=4$ in
this paper. Furthermore, each element will be associated with 3 additional
parameters: the thickness parameter ${\beta}^{(e)}$, the anisotropy parameter
${\alpha}^{(e)}$, and the orientation parameter ${\theta}^{(e)}$. Observe that
the parameters associated with an element are common to all sites within the
element. By optimizing all parameters over all elements (as discussed in the
remainder of the paper), one can design a multi-scale porous structure; this
is schematically illustrated in fig. 6. However, direct optimization will
still be computationally prohibitive due to the need for high resolution
finite element analysis.

Figure 6: Given a discretized domain, we wish to populate each element with
Voronoi cells. All cells in a particular microstructure have the same
thickness ($\beta$), and are anisotropic ($\alpha$) with an orientation
($\theta$).

To reduce the computational cost we propose a two-stage strategy (see the
graphical abstract fig. 1):

  1. 1.

Offline Computation: This involves generating numerous representative Voronoi
microstructures, computing their homogenized properties and training of NNs
(see section 2.4).

  2. 2.

Multi-Scale Optimization: Exploiting the trained NNs to carry out efficient
multi-scale Voronoi structure generation (see section 2.5).

###  2.4 Offline Computation

We now describe various steps in the offline computation.

####  2.4.1 Voronoi Microstructure Generation

The first step is to create a large, representative set of Voronoi
microstructures associated with the macro elements. A naive approach is to
generate $S$ random sites (points) within a macro element and generate the
corresponding Voronoi microstructure. However, this would not represent a
typical microstructure since points in the neighboring macro elements also
influence the Voronoi microstructure. For example, the naive approach would
fail to generate solid material on the element’s boundary. We therefore
consider a macro element and its 8 neighboring macro elements, and generate
random cell sites in all 9 elements; see fig. 7a. Then, a Voronoi
microstructure is generated in the central element using all the cell sites,
and a random set of parameters [67] $\beta$, $\alpha$ and $\theta$, that are
generated uniformly over a pre-defined range, as described in numerical
experiments; see fig. 7b for a typical Voronoi microstructure.

Figure 7: Data generation involves (a) distributing Voronoi cell site
locations for the central microstructure element (1) and its neighbors (8),
(b) generating Voronoi microstructures based on cell locations, thickness,
degree of anisotropy and rotation values for the central microstructure
element (c) extracting the central microstructure design.

Further, to ensure that the points are sufficiently separated from one
another, we model the points as a perturbation $(\Delta_{x},\Delta_{y})$ from
grid points; see fig. 8.

Figure 8: The cell sites
$(\mathcolor{red}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}})$
within a microstructure are obtained as a perturbation to a grid of points
$(\boldsymbol{\bigstar})$.

####  2.4.2 Homogenization

Once these Voronoi microstructures are generated, the macro element is
discretized into a grid of $120\times 120$ _micro_ elements; see fig. 9a. The
mesh is then used for numerical homogenization to compute the elasticity
matrix $\mathbf{C}$ and the volume fraction $v$; see fig. 9b and fig. 9c.
Without a loss of generality, we assume a Young’s Modulus $E=1$ and a
Poisson’s ratio $\nu=0.3$; see [68] for details on numerical homogenization.

Figure 9: Numerical homogenization:(a) Discretization of the density field
onto a mesh, (d) employing numerical homogenization to determine the (c)
constitutive properties of the microstructure.

####  2.4.3 Cholesky Decomposition

In the subsequent sections, we employ a NN to establish a mapping between the
Voronoi parameters, and the corresponding homogenized elasticity matrix
components (and volume fraction). However, a direct mapping might lead to the
prediction of stiffness matrices that are not positive definite [69]. To avoid
this, we first carry out a Cholesky decomposition of $\mathbf{C}$:

$\mathbf{C}=\mathbf{L}\mathbf{L}^{T}$ (6)

where,

$\mathbf{L}=\begin{bmatrix}L_{00}&0&0\\\ L_{10}&L_{11}&0\\\
L_{20}&L_{21}&L_{22}\\\ \end{bmatrix}$

The components of $\mathbf{L}$, i.e., $\\{L_{00},L_{10}...,L_{22}\\}$ are then
used for training the NN, as discussed next.

####  2.4.4 Neural Network Training

Following the data generation and Cholesky decomposition, we use NNs [70] to
establish a mapping between the Voronoi parameters and computed data. NNs have
been utilized for various material modeling applications, including modeling
the plastic response of metal [71], to construct a relation between
macroscopic stress and crack opening responses [72], homogenization of
composite structures [73, 74], constitutive modeling of elastomeric foams
[75], and nonlinear response of carbon nanotubes [76]. In addition, NNs have
been utilized in the context of multi-scale topology optimization [77].

For the design of microstructures using NNs pixel-based representations have
been employed [78, 79, 80]. This method assumes designs can be constructed
using solid or void elements in space. While this allows extensive design
flexibility, it has drawbacks [81]. There is a considerable computational
burden for both design evaluation and machine learning due to limitations in
scalability with respect to resolution [81]. Second, efficient exploration of
the design space requires the utilization of constraints to achieve desired
design attributes [82]. Many studies [83, 84, 85, 86] employing NNs trained on
pixel-based representations to predict microstructure properties typically
optimize the networks to minimize the difference between the predicted design
and a ground truth solution using loss functions such as mean squared error
(MSE), mean absolute error (MAE), or binary cross-entropy. However, this
approach may result in inaccurate predictions as structurally similar designs
in a pixelated format can exhibit significantly different properties [81, 87].

Here, NNs are employed to predict the homogenized elastic response (and volume
fraction) as a continuous and differentiable mapping of microstructure
parameters, facilitating gradient-based optimization [88]. The proposed NN
architecture consists of the following components (see fig. 10):

  1. 1.

The input consists of the $(x,y)$ location of the sites (points) from all 9
macro elements (the elements are ordered bottom to top, then left to right,
and the parameters $\beta$, $\alpha$ and $\theta$, associated with the center
element, i.e., the input is a ($9\times 2\times S+3$) dimensional vector,
where $S$ is the number of cell sites per element.

  2. 2.

The NN is a fully connected network consisting of four hidden layers, each
comprising $50$ neurons with ReLU activation functions.

  3. 3.

The output layer consists the $6$ components of $\mathbf{\hat{L}}$ (the
predicted values of $\mathbf{L}$) and the volume fraction $\hat{v}$ (the
predicted value of $v$), i.e., it is a $7$ dimensional vector denoted as
$\mathbf{\hat{\Psi}}$ (the predicted value of $\mathbf{\Psi}$). To ensure that
the element stiffness matrix remains positive definite, the predicted diagonal
components $\\{\hat{L}^{(e)}_{00},\hat{L}^{(e)}_{11},\hat{L}^{(e)}_{22}\\}$
are clipped, from below, to a small positive value (of $10^{-6}$).

Figure 10: Proposed surrogate neural network architecture.

The NN is then trained on the $(N_{d})$ generated data, i.e., the NN minimizes
the mean squared difference between the actual and predicted microstructure
properties. The loss function is defined here to be:

$L_{NN}=\frac{1}{N_{d}}\big{(}||\mathbf{\hat{\Psi}}-\mathbf{\Psi}||^{2}_{2}\big{)}$
(7)

###  2.5 Multi-Scale Optimization

We now consider the proposed multi-scale optimization framework illustrated in
fig. 11. Recall that each macro element is associated with $S$ cell cites and
3 geometric parameters $\beta^{(e)}$, $\theta^{(e)}$ and $\alpha^{(e)}$. These
cell sites and geometric parameters will be optimized to meet a desired
objective. To avoid expensive homogenization, the trained NN can be exploited
to predict the elasticity response of each element. Specifically, at each step
of the optimization process, the current locations of the cell sites of $e$
and its 8 neighbors, together with the current values of 3 geometric
parameters of $e$ are extracted. We then apply a smoothing radial filter on
the three parameters, i.e., at each element, a weighted average of every
parameter about the neighborhood of that element is computed. This filtering
process helps promote connectivity by smoothing out abrupt changes across
neighboring microstructures.

Using these smoothed parameters, the trained NN is used to estimate
$\mathbf{\hat{L}}^{(e)}$, and therefore $\mathbf{\hat{C}}^{(e)}$. Given the
elasticity matrix $\mathbf{\hat{C}}^{(e)}$ of a macro element, recall that the
element stiffness matrix is defined as [89]:

$[K^{(e)}]=\int\limits_{\Omega^{(e)}}[\nabla
N^{(e)}]^{T}[\mathbf{\hat{C}}^{(e)}][\nabla N^{(e)}]d\Omega^{(e)}$ (8)

where $[\nabla N^{(e)}]$ is the gradient of the shape matrix. Once the element
stiffness matrices are computed, the global stiffness matrix $\bm{K}$ and
force vector $\bm{f}$ are assembled, followed by the solution of the
displacement field $\bm{u}=\bm{K}^{-1}\bm{f}$. The objective, etc, are then
computed followed by an update of all design parameters. These steps are
discussed in detail in the remainder of this section.

Figure 11: The multiscale topology optimization framework utilizes the cell
site locations of the microstructure and its neighboring microstructures,
along with the microstructure’s thickness, degree of anisotropy, and
orientation to predict the constitutive properties. This prediction is
achieved using the trained neural networks, facilitating finite element
analysis in multiscale optimization. Note: Each microstructure element
comprises $S$ sites, and for a neighborhood of 9 elements, the neural network
is provided with an input of $9S$ sites.

####  2.5.1 Optimization Problem

We now describe the multi-scale optimization formulation.

Design Variables: The design variables for each element are represented by
$\bm{\zeta}^{(e)}=\\{\bm{\Delta}_{x}^{(e)},\bm{\Delta}_{y}^{(e)},\beta^{(e)},\alpha^{(e)},\theta^{(e)}\\}$.
The complete set of design variables is represented as
$\bm{\overline{\zeta}}=\\{\bm{\zeta}^{(1)},\bm{\zeta}^{(2)},....,\bm{\zeta}^{(n)}\\}$.
where $n$ is the number of macro elements.

Objective: We consider a simple compliance minimization objective:

$J(\overline{\bm{\zeta}})=\bm{f}^{T}\bm{u}$ (9)

where $\bm{u}$ the nodal displacements, and $\bm{f}$ is the imposed load.

Volume Constraint: A global volume constraint is imposed:

$g_{V}(\overline{\bm{\zeta}})\equiv\frac{1}{n}\frac{\sum\limits_{e=1}^{n}\hat{v}^{(e)}}{v_{max}}-1\leq
0$ (10)

where $\hat{v}^{(e)}$ is the predicted volume fraction, and $v_{max}$ is the
imposed upper bound on the global volume fraction.

Consequently, the optimization problem can be formulated as follows:

$\displaystyle\underset{\overline{\bm{\zeta}}=\\{\bm{\zeta}_{1},\bm{\zeta}_{2},\ldots\bm{\zeta}_{n}\\}}{\text{minimize}}$
$\displaystyle J(\overline{\bm{\zeta}})=\bm{f}^{T}\bm{u}$ (11a) subject to
$\displaystyle\bm{K}(\overline{\bm{\zeta}})\bm{u}=\bm{f}$ (11b) and
$\displaystyle
g_{V}(\overline{\bm{\zeta}})\equiv\frac{\sum\limits_{e=1}^{n}\hat{v}^{(e)}}{nv_{max}}-1\leq
0$ (11c)

####  2.5.2 Loss Function

The constrained minimization problem eq. 11 is transformed into an
unconstrained loss function minimization, using the penalty scheme [90].
Specifically, the loss function is defined as:

$\mathcal{L}(\bm{\zeta})=\frac{J(\bm{\zeta})}{J^{0}}+\gamma
g_{v}(\bm{\zeta})^{2}$ (12)

where $J^{0}$ is the initial compliance. The constraint penalty parameter
$\gamma$ begins with an initial value of $\gamma=0.1$ and is subsequently
incremented by $\Delta\gamma=0.25$ after each epoch. The gradient-based Adam
optimizer [91] is used to minimize eq. 12.

####  2.5.3 Sensitivity

A critical aspect of gradient-based optimization is determining the
sensitivity, or derivatives, of both the objective function and constraints
with respect to the optimization parameters. Here, we harness the automatic
differentiation (AD) capabilities of the PyTorch framework [92] to avoid
manual sensitivity calculations [93]. In practical terms, we only need to
define the forward expressions, and PyTorch’s computational library will
compute all necessary derivatives with machine precision [94].

##  3 Numerical Experiments

In this section, we conduct several experiments to illustrate the method.
Without loss of generality, the default parameters are as follows:

  1. 1.

We assume the domain to consist of $40\times 20$ elements, where each element
is assumed to be of unit area.

  2. 2.

The material is assumed to have Young’s modulus of 1 and Poisson’s ratio of
0.3.

  3. 3.

The force is assumed to be 1 unit.

  4. 4.

Termination criteria include an optimization loss tolerance of $10^{-3}$ or a
maximum of 300 iterations.

  5. 5.

All experiments were conducted on a MacBook M2 Air with $16$ GB RAM.

Other parameters are provided in the corresponding section.

###  3.1 Offline Experiments

####  3.1.1 Dataset Generation

To establish a mapping between microstructure element parameters and
microstructure properties using a surrogate NN, the first step involves
acquiring training data. This process entails allowing the $4$ seed points per
microstructure element $(\bar{x},\bar{y})$ to be uniformly distributed within
a unit-length element. This distribution is achieved by varying the
perturbation from the base grid points $(\Delta_{x},\Delta_{y})$ within the
range of $[-0.225,0.225]$, ensuring a minimum separation of $0.1$ between
neighboring cell sites.

Additionally, the parameters thickness ($\beta$), anisotropy ($\alpha$) and
orientation ($\theta$) are considered to be in range $[0.3,3]$, $[1,3.5]$ and
$[0,\pi]$ respectively. These ranges are chosen to encompass a wide variation
of volume fractions observed across $12000$ samples of microstructure elements
(as depicted in fig. 12 (a)) and to account for anisotropic homogenized
elasticity matrix (as shown in fig. 12 (b)). The thickness ($\beta$),
anisotropy ($\alpha$) and orientation ($\theta$) were filtered using a filter
radius of 3.

Figure 12: Distribution of (a) volume fraction and (b) constitutive matrix
components $C_{00}$ and $C_{11}$ as a function of microstructure density.

####  3.1.2 NN Training

The training data generated in the previous section serves as input for the
surrogate NN, which consists of an architecture outlined in section 2.4.4. The
network’s output comprises the microstructure volume fraction and Cholesky
factors of the homogenized elasticity matrix. The data is split into three
sets: training ($N_{d}=10000$ data points), validation ($1000$ data points),
and testing ($1000$ data points). During training, the NN aims to minimize the
mean squared error between its predicted and actual microstructure properties.
This training process is performed for a maximum of $300$ iterations, with a
learning rate of $5*10^{-5}$ and a batch size of $64$. The training process is
stopped at the maximum iteration number or when the validation loss starts to
increase, thus preventing over-fitting. Following this approach, we observe
convergence in training, test, and validation loss between actual and
predicted volume fraction and homogenized elasticity matrix, as depicted in
fig. 13.

Figure 13: Convergence of the training, testing, and validation loss.

###  3.2 Multi-Scale Optimization Experiments

Next, we demonstrate multi-scale optimization through several experiments.

####  3.2.1 Validation: Tensile Bar

Consider the tensile bar problem in fig. 14(a). The objective is to find the
optimal topology that minimizes compliance subjected to a volume fraction
$v=0.4$ as described in [89]. The compliance of single-scale optimized design
reported in reference [89] is $183$. Here, we obtain a multiscale optimized
design, illustrated in fig. 15 (c). As one can observe the porous structure
resembles an ideal tensile bar, while meeting the minimum porosity imposed
(via the parameter $\beta$). Furthermore, near the transition, the porous
structures are oriented along the axis, as expected. The compliance of the
porous structure is $209$.

Figure 14: Validation of the multiscale framework: (a) tensile bar boundary
conditions, (b) multiscale design [89] and (c) multiscale porous design.

####  3.2.2 Comparison

Next consider the mid-cantilever problem in fig. 15(a). The objective is to
find the optimal topology that minimizes compliance subject to a volume
fraction $v=0.6$. The single-scale optimized design achieved using the code
from [95] is illustrated in fig. 15(b); the final compliance is $61.05$. In
[96], by imposing a local volume constraint of $v=0.6$, an optimized porous
structure with compliance of $76.86$ was reported. Here, we obtain a
multiscale optimized design, illustrated in fig. 15 (c), with a compliance of
$68.3$. The computation took $13.7$ seconds. Once again, the computed
structure resembles the single-scale structure but deviates from it to meet
the porosity constraints.

To evaluate the accuracy of the NN mapping, we reconstructed the
microstructures of the optimized design. We then computed the actual
homogenized matrices via FEA of each macro element. Subsequent global analysis
yielded a compliance of $74.6$ and a volume fraction of $0.65$. This
translates to an error of $8.4\%$ and $6.1\%$ for the compliance and volume
fraction, respectively.

Figure 15: Validation of the multiscale framework: (a) mid cantilever boundary
conditions, (b) single scale design and (c) multiscale porous design.

####  3.2.3 Convergence against Single-Scale

We illustrate the typical convergence of the proposed algorithm for the mid-
cantilever problem in fig. 15(a). The compliance and volume constraint is
illustrated in fig. 16. Similar convergence behavior was observed for all
other examples.

Figure 16:  Convergence of the objective and constraint during optimization.

####  3.2.4 Parameter Variations

A central hypothesis of our current work is that better multiscale designs can
be obtained with a broader range of parameters: thickness, anisotropy, and
orientation. To validate these hypotheses, we revisit the problem depicted in
fig. 15(a), enforcing a volume constraint of $0.5$ but considering various
restrictions on the parameters.

Thickness Parameter: $\beta$

In this numerical experiment, we set the lower bound of the thickness
parameter at $0.3$ and analyze the impact of thickness on both the objective
function and computational time by setting its upper bound to $1.$, $2.$, and
$3.$. The resulting topologies, illustrated in fig. 17, confirm our
expectations: the objective function improves with an increasing range of
thickness parameters. Furthermore, the computational time was approximately
$13.7s$, regardless of the range of thickness parameter.

Figure 17: Impact of thickness parameter on compliance.

Anisotropy Parameter: $\alpha$

In this numerical experiment, we set the lower bound for the degree of
anisotropy at $1$ and investigate the influence of anisotropy on both the
objective function and computational time by setting its upper bound to $1.5$,
$2.5$, and $3.5$. The resulting topologies, as depicted in fig. 18, align with
our anticipated outcomes: the objective decreases with increasing maximum
degrees of anisotropy. Once again, the computational time was unaffected by
the change in upper bound.

Figure 18: Varying maximum degree of anisotropy with minimum degree of
anisotropy = $1$

Orientation Parameter: $\theta$

Next, we kept the orientation fixed at $\theta=0$ and optimized the design. As
expected, our observations, fixing the orientation results in lower
performance as seen in fig. 19.

Figure 19: Impact of orientation on compliance.

Volume Fraction Variation

Finally, we consider the trade-off between the objective (compliance) and
constraint (volume fraction) through exploration of the Pareto front is
crucial for making informed design choices. In this study, we examined the
heel bone problem depicted in fig. 3 (a). We determined the optimal topologies
for various volume fraction constraints. The results are illustrated in fig.
20; as expected, the compliance increases as the volume fraction decreases. To
evaluate the accuracy of the NN mapping, we recomputed the actual values using
a macroscale FEA/ homogenization approach. The error in compliance and volume
fraction was less than $10\%$ in all cases.

Figure 20: Trade-off between compliance and volume fraction.

###  3.3 Is Offline Computation Justified?

The offline homogenization and data generation of $12000$ microstructures
required $26.4$ minutes, while the training of the surrogate NN required $2$
minutes. Thus, the total one-time overhead amounted to $28.4$ minutes.
Subsequently, the optimization of the mid-cantilever, for example, (using a
grid size of $40\times 20$) consumed $0.24$ minutes.

Now, let’s consider two hypothetical scenarios. The first scenario is a brute-
force fine-scale optimization. We note that the grid size of $40\times 40$
microscale elements is required for each homogenization of the microstructure
element. For the single-scale optimization of the mid-cantilever (using a grid
size of $40\times 20$ of microstructure elements), one requires a grid of
$40\times 40\times 40\times 20$, which is $1600\times 800$ elements. One
iteration of single-scale optimization with $1600\times 800$ elements takes
$16.62$ minutes. Consequently, the anticipated total optimization time amounts
to $16.62\times 300$, equating to $4986$ minutes.

Next, consider a scenario where we carry out multiscale optimization but do
not rely on offline training, i.e., we carry out concurrent homogenization-
based multiscale optimization. Observe that the time required for each
homogenization is $26.4/12000$, i.e., $0.132$ seconds. Now, for concurrent
homogenization of the mid-cantilever, the task entails homogenization across
every element within a $40\times 20$ grid over $300$ iterations. As a result,
the projected optimization duration can be estimated as $0.132\times 800\times
300/60$, which equals $528$ minutes.

Thus, the proposed offline NN-based multiscale optimization is computationally
far superior. Furthermore, the loss in accuracy in predicted compliance and
volume fraction was found to be less than $10\%$ in all experiments.

##  4 Conclusion

In this paper, we introduced a novel multi-scale topology optimization
framework using Voronoi microstructures. Our approach involves an offline
homogenization process coupled with NN training to establish a continuous and
differentiable mapping of microstructural geometric parameters to constitutive
properties. Subsequently, we carried out multiscale optimization to minimize
compliance subject to volume constraints.

Our numerical results illustrate that the proposed method offers significant
computational advantages over concurrent homogenization, with minimal loss in
accuracy. Additionally, we observe that Voronoi microstructures, along with
parameters such as thickness anisotropy and orientation, expand the design
space. Furthermore, considering neighboring microstructure cell sites during
NN training facilitates macroscopic connectivity. Finally, training the NN
with Cholesky factors ensures the positive definiteness of the constitutive
matrix.

Future research avenues include extending the framework to encompass Voronoi
structures found in natural systems, which often serve multiple functions. For
example, in bone infills, pores act as conduits for transport while the
material provides structural integrity. This could involve a multi-objective
formulation integrating compliance and diffusivity to generate bone-like
porous structures. In natural porous materials, the density of Voronoi cell
sites varies, whereas in the present formulation, the number of cell sites in
a region remains constant. Extending the framework to three dimensions is also
a crucial area for future investigation. Additionally, the proposed framework
focused on a singular class of microstructures derived from Turing patterns,
but the potential for generating a broader spectrum of patterns exists through
the application of reaction-diffusion equations [97].

## Acknowledgments

The University of Wisconsin, Madison Graduate School supported this work.

## Compliance with ethical standards

The authors declare that they have no conflict of interest.

## Replication of Results

The Python code is available at github.com/UW-ERSL/VoroTO

## References

  * [1] Martin Philip Bendsoe and Ole Sigmund.  Topology optimization: theory, methods, and applications.  Springer Science & Business Media, 2013.
  * [2] Ole Sigmund and Kurt Maute.  Topology optimization approaches: A comparative review.  Structural and multidisciplinary optimization, 48(6):1031–1055, 2013.
  * [3] Ronald F Gibson.  A review of recent research on mechanics of multifunctional composite materials and structures.  Composite structures, 92(12):2793–2810, 2010.
  * [4] Franz Knoll and Thomas Vogel.  Design for robustness, volume 11.  IABSE, 2009.
  * [5] Zhong Hu, Kaushik Thiyagarajan, Amrit Bhusal, Todd Letcher, Qi Hua Fan, Qiang Liu, and David Salem.  Design of ultra-lightweight and high-strength cellular structural composites inspired by biomimetics.  Composites Part B: Engineering, 121:108–121, 2017.
  * [6] M Sadegh Ebrahimi, R Hashemi, and E Etemadi.  In-plane energy absorption characteristics and mechanical properties of 3d printed novel hybrid cellular structures.  Journal of Materials Research and Technology, 20:3616–3632, 2022.
  * [7] Luis Cardoso, Susannah P Fritton, Gaffar Gailani, Mohammed Benalla, and Stephen C Cowin.  Advances in assessment of bone porosity, permeability and interstitial fluid flow.  Journal of biomechanics, 46(2):253–265, 2013.
  * [8] Michael F Ashby and Lorna J Gibson.  Cellular solids: structure and properties.  Press Syndicate of the University of Cambridge, Cambridge, UK, pages 175–231, 1997.
  * [9] Helge Klippstein, Hany Hassanin, Alejandro Diaz De Cerio Sanchez, Yahya Zweiri, and Lakmal Seneviratne.  Additive manufacturing of porous structures for unmanned aerial vehicles applications.  Advanced Engineering Materials, 20(9):1800290, 2018.
  * [10] Shiqi Wu, Daming Chen, Guangdong Zhao, Yuan Cheng, Boqian Sun, Xiaojie Yan, Wenbo Han, Guiqing Chen, and Xinghong Zhang.  Controllable synthesis of a robust sucrose-derived bio-carbon foam with 3d hierarchical porous structure for thermal insulation, flame retardancy and oil absorption.  Chemical Engineering Journal, 434:134514, 2022.
  * [11] Scott J Hollister.  Porous scaffold design for tissue engineering.  Nature materials, 4(7):518–524, 2005.
  * [12] Young-Joon Seol, Dong Yong Park, Ju Young Park, Sung Won Kim, Seong Jin Park, and Dong-Woo Cho.  A new method of fabricating robust freeform 3d ceramic scaffolds for bone tissue regeneration.  Biotechnology and Bioengineering, 110(5):1444–1455, 2013.
  * [13] Da Chen, Kang Gao, Jie Yang, and Lihai Zhang.  Functionally graded porous structures: Analyses, performances, and applications–a review.  Thin-Walled Structures, 191:111046, 2023.
  * [14] Diplesh Gautam and Venkatesh KP Rao.  Nondestructive evaluation of mechanical properties of femur bone.  Journal of Nondestructive Evaluation, 40(1):22, 2021.
  * [15] Chukwuzubelu Okenwa Ufodike, Mohammad Faisal Ahmed, and Grzegorz Dolzyk.  Additively manufactured biomorphic cellular structures inspired by wood microstructure.  Journal of the Mechanical Behavior of Biomedical Materials, 123:104729, 2021.
  * [16] SR Jongerius and D Lentink.  Structural analysis of a dragonfly wing.  Experimental Mechanics, 50:1323–1334, 2010.
  * [17] Haizum Aimi Zaharin, Ahmad Majdi Abdul Rani, Farooq I Azam, Turnad Lenggo Ginta, Nabihah Sallih, Azlan Ahmad, Nurul Azhani Yunus, and Tun Zainal Azni Zulkifli.  Effect of unit cell type and pore size on porosity and mechanical behavior of additively manufactured ti6al4v scaffolds.  Materials, 11(12):2402, 2018.
  * [18] Christopher R Jacobs, Juan C Simo, Gary S Beaupre, and Dennis R Carter.  Adaptive bone remodeling incorporating simultaneous density and anisotropy considerations.  Journal of biomechanics, 30(6):603–613, 1997.
  * [19] Jaejong Park, Alok Sutradhar, Jami J Shah, and Glaucio H Paulino.  Design of complex bone internal structure using topology optimization with perimeter control.  Computers in biology and medicine, 94:74–84, 2018.
  * [20] Yohann Bala, Emmanuelle Lefèvre, Jean-Paul Roux, Cécile Baron, Philippe Lasaygues, Martine Pithioux, Valérie Kaftandjian, and Hélène Follet.  Pore network microarchitecture influences human cortical bone elasticity during growth and aging.  Journal of the mechanical behavior of biomedical materials, 63:164–173, 2016.
  * [21] Ibrahim Goda, Jean-Francois Ganghoffer, Sławomir Czarnecki, Radosław Czubacki, and Paweł Wawruch.  Topology optimization of bone using cubic material design and evolutionary methods based on internal remodeling.  Mechanics Research Communications, 95:52–60, 2019.
  * [22] Phong Cong Hong Nguyen, Youngdoo Kim, and Young Choi.  Lightweight design with metallic additively manufactured cellular structures.  Journal of Computational Design and Engineering, 9(1):155–167, 2022.
  * [23] Wahyudin P Syam, Wu Jianwei, Bo Zhao, Ian Maskery, Waiel Elmadih, and Richard Leach.  Design and analysis of strut-based lattice structures for vibration isolation.  Precision Engineering, 52:494–506, 2018.
  * [24] Lorenzo Valdevit, Alan J Jacobsen, Julia R Greer, and William B Carter.  Protocols for the optimal design of multi-functional cellular structures: from hypersonics to micro-architected materials.  Journal of the American Ceramic Society, 94:s15–s34, 2011.
  * [25] Gang Lei, PC Dong, ZS Wu, SH Gai, SY Mo, and Z Li.  Multi-scale structures of porous media and the flow prediction.  Journal of Natural Gas Science and Engineering, 21:986–992, 2014.
  * [26] Jonàs Martínez, Jérémie Dumas, and Sylvain Lefebvre.  Procedural voronoi foams for additive manufacturing.  ACM Transactions on Graphics (TOG), 35(4):1–12, 2016.
  * [27] Guanjun Wang, Lida Shen, Jianfeng Zhao, Huixin Liang, Deqiao Xie, Zongjun Tian, and Changjiang Wang.  Design and compressive behavior of controllable irregular porous scaffolds: based on voronoi-tessellation and for additive manufacturing.  ACS biomaterials science & engineering, 4(2):719–727, 2018.
  * [28] Clément Audibert, Julien Chaves-Jacob, Jean-Marc Linares, and Quentin-Alexis Lopez.  Bio-inspired method based on bone architecture to optimize the structure of mechanical workspieces.  Materials & Design, 160:708–717, 2018.
  * [29] Jianming Ying, Lin Lu, Lihao Tian, Xin Yan, and Baoquan Chen.  Anisotropic porous structure modeling for 3d printed objects.  Computers & Graphics, 70:157–164, 2018.
  * [30] Jun Wu, Ole Sigmund, and Jeroen P Groen.  Topology optimization of multi-scale structures: a review.  Structural and Multidisciplinary Optimization, 63:1455–1480, 2021.
  * [31] Rahul Kumar Padhy, Aaditya Chandrasekhar, and Krishnan Suresh.  Fluto: Graded multi-scale topology optimization of large contact area fluid-flow devices using neural networks.  Engineering with Computers, pages 1–17, 2023.
  * [32] Rahul Kumar Padhy, Krishnan Suresh, and Aaditya Chandrasekhar.  Tomas: Topology optimization of multiscale fluid devices using variational autoencoders and super-shapes.  arXiv preprint arXiv:2309.08435, 2023.
  * [33] Jun Wu, Christian Dick, and Rüdiger Westermann.  A system for high-resolution topology optimization.  IEEE transactions on visualization and computer graphics, 22(3):1195–1208, 2015.
  * [34] Jun Wu, Weiming Wang, and Xifeng Gao.  Design and optimization of conforming lattice structures.  IEEE transactions on visualization and computer graphics, 27(1):43–56, 2019.
  * [35] James K Guest.  Imposing maximum length scale in topology optimization.  Structural and Multidisciplinary Optimization, 37:463–473, 2009.
  * [36] Suguang Dou.  A projection approach for topology optimization of porous structures through implicit local volume control.  Structural and Multidisciplinary Optimization, 62(2):835–850, 2020.
  * [37] Sourav Das and Alok Sutradhar.  Multi-physics topology optimization of functionally graded controllable porous structures: Application to heat dissipating problems.  Materials & Design, 193:108775, 2020.
  * [38] Zhi Zhao and Xiaojia Shelly Zhang.  Design of graded porous bone-like structures via a multi-material topology optimization approach.  Structural and Multidisciplinary Optimization, 64:677–698, 2021.
  * [39] Fan Feng, Shiying Xiong, Ziyue Liu, Zangyueyang Xian, Yuqing Zhou, Hiroki Kobayashi, Atsushi Kawamoto, Tsuyoshi Nomura, and Bo Zhu.  Cellular topology optimization on differentiable voronoi diagrams.  International Journal for Numerical Methods in Engineering, 124(1):282–304, 2023.
  * [40] Hongjia Lu, Ting-Uei Lee, Jiaming Ma, Da Chen, and Yi Min Xie.  Designing 2d stochastic porous structures using topology optimisation.  Composite Structures, 321:117305, 2023.
  * [41] Yanfang Zhao, Guikai Guo, Xinyu Xie, and Wenjie Zuo.  Concurrent multiscale topology optimization of hollow structures considering geometrical nonlinearity.  Engineering with Computers, pages 1–18, 2023.
  * [42] Martin P Bendsøe and Ole Sigmund.  Material interpolation schemes in topology optimization.  Archive of applied mechanics, 69:635–654, 1999.
  * [43] Jun Yan, Qianqian Sui, Zhirui Fan, Zunyi Duan, and Tao Yu.  Clustering-based multiscale topology optimization of thermo-elastic lattice structures.  Computational Mechanics, 66:979–1002, 2020.
  * [44] Hao Li, Zhen Luo, Nong Zhang, Liang Gao, and Terry Brown.  Integrated design of cellular composites using a level-set topology optimization method.  Computer Methods in Applied Mechanics and Engineering, 309:453–475, 2016.
  * [45] Christian Schumacher, Steve Marschner, Markus Gross, and Bernhard Thomaszewski.  Mechanical characterization of structured sheet materials.  ACM Transactions on Graphics (TOG), 37(4):1–15, 2018.
  * [46] Davi Colli Tozoni, Jérémie Dumas, Zhongshi Jiang, Julian Panetta, Daniele Panozzo, and Denis Zorin.  A low-parametric rhombic microstructure family for irregular lattices.  ACM Transactions on Graphics (TOG), 39(4):101–1, 2020.
  * [47] Xin Yan, Cong Rao, Lin Lu, Andrei Sharf, Haisen Zhao, and Baoquan Chen.  Strong 3d printing by tpms injection.  IEEE Transactions on Visualization and Computer Graphics, 26(10):3037–3050, 2019.
  * [48] Jiangbei Hu, Shengfa Wang, Baojun Li, Fengqi Li, Zhongxuan Luo, and Ligang Liu.  Efficient representation and optimization for tpms-based porous structures.  IEEE Transactions on Visualization and Computer Graphics, 28(7):2615–2627, 2020.
  * [49] Ming Li, Jingqiao Hu, Wei Chen, Weipeng Kong, and Jin Huang.  Explicit topology optimization of conforming voronoi foams.  arXiv preprint arXiv:2308.04001, 2023.
  * [50] Bruno Jetté, Vladimir Brailovski, Mathieu Dumas, Charles Simoneau, and Patrick Terriault.  Femoral stem incorporating a diamond cubic lattice structure: Design, manufacture and testing.  Journal of the Mechanical Behavior of Biomedical Materials, 77:58–72, 2018.
  * [51] Farnaz Ghorbani, Dejian Li, Shuo Ni, Ying Zhou, and Baoqing Yu.  3d printing of acellular scaffolds for bone defect regeneration: A review.  Materials Today Communications, 22:100979, 2020.
  * [52] Long Chao, Yangdong He, Jiasen Gu, Deqiao Xie, Youwen Yang, Lida Shen, Guofeng Wu, Lin Wang, and Zongjun Tian.  Evaluation of compressive and permeability behaviors of trabecular-like porous structure with mixed porosity based on mechanical topology.  Journal of Functional Biomaterials, 14(1):28, 2023.
  * [53] Hong-Yuan Lei, Jing-Rong Li, Zhi-Jia Xu, and Qing-Hui Wang.  Parametric design of voronoi-based lattice porous structures.  Materials & Design, 191:108607, 2020.
  * [54] Lin Lu, Andrei Sharf, Haisen Zhao, Yuan Wei, Qingnan Fan, Xuelin Chen, Yann Savoye, Changhe Tu, Daniel Cohen-Or, and Baoquan Chen.  Build-to-last: Strength to weight 3d printed objects.  ACM Transactions on Graphics (ToG), 33(4):1–10, 2014.
  * [55] XY Kou and ST Tan.  A simple and effective geometric representation for irregular porous structure modeling.  Computer-Aided Design, 42(10):930–941, 2010.
  * [56] Zhongren Wang, Lihao Tian, Xiaokang Liu, Andrei Sharf, and Lin Lu.  Stochastic porous microstructures.  arXiv preprint arXiv:2305.09176, 2023.
  * [57] Filippo Cucinotta, Marcello Raffaele, and Fabio Salmeri.  A stress-based topology optimization method by a voronoi tessellation additive manufacturing oriented.  The International Journal of Advanced Manufacturing Technology, 103:1965–1975, 2019.
  * [58] Lianxiong Chen, Ran Zhang, Xihua Chu, and Hui Liu.  Multiscale design of graded stochastic cellular structures for the heat transfer problem.  Applied Sciences, 13(7):4409, 2023.
  * [59] Der-Tsai Lee and Robert L Drysdale, III.  Generalization of voronoi diagrams in the plane.  SIAM Journal on Computing, 10(1):73–87, 1981.
  * [60] Howie Choset and Joel Burdick.  Sensor-based exploration: The hierarchical generalized voronoi graph.  The International Journal of Robotics Research, 19(2):96–125, 2000.
  * [61] Kenneth E Hoff III, John Keyser, Ming Lin, Dinesh Manocha, and Tim Culver.  Fast computation of generalized voronoi diagrams using graphics hardware.  In Proceedings of the 26th annual conference on Computer graphics and interactive techniques, pages 277–286, 1999.
  * [62] Elena Arseneva and Evanthia Papadopoulou.  Randomized incremental construction for the hausdorff voronoi diagram revisited and extended.  Journal of Combinatorial Optimization, 37:579–600, 2019.
  * [63] Maxime Sainlot, Vincent Nivoliers, and Dominique Attali.  Restricting voronoi diagrams to meshes using corner validation.  In Computer Graphics Forum, volume 36, pages 81–91. Wiley Online Library, 2017.
  * [64] Francois Labelle and Jonathan Richard Shewchuk.  Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation.  In Proceedings of the nineteenth annual symposium on Computational geometry, pages 191–200, 2003.
  * [65] TFW Van Nuland, JAW Van Dommelen, and MGD Geers.  An anisotropic voronoi algorithm for generating polycrystalline microstructures with preferred growth directions.  Computational Materials Science, 186:109947, 2021.
  * [66] Azwirman Gusrialdi, Sandra Hirche, Takeshi Hatanaka, and Masayuki Fujita.  Voronoi based coverage control with anisotropic sensors.  In 2008 American control conference, pages 736–741. IEEE, 2008.
  * [67] Travis E Oliphant et al.  Guide to numpy, volume 1.  Trelgol Publishing USA, 2006.
  * [68] Erik Andreassen and Casper Schousboe Andreasen.  How to determine composite material properties using numerical homogenization.  Computational Materials Science, 83:488–495, 2014.
  * [69] Charles F Jekel, Kenneth E Swartz, Daniel A White, Daniel A Tortorelli, and Seth E Watts.  Neural network layers for prediction of positive definite elastic stiffness tensors.  arXiv preprint arXiv:2203.13938, 2022.
  * [70] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366, 1989.
  * [71] Tomonari Furukawa and Mark Hoffman.  Accurate cyclic plastic analysis using a neural network material model.  Engineering Analysis with Boundary Elements, 28(3):195–204, 2004.
  * [72] Jörg F Unger and Carsten Könke.  Neural networks as material models within a multiscale approach.  Computers & structures, 87(19-20):1177–1186, 2009.
  * [73] M Lefik, DP Boso, and BA Schrefler.  Artificial neural networks in numerical modelling of composites.  Computer Methods in Applied Mechanics and Engineering, 198(21-26):1785–1804, 2009.
  * [74] Hou Man and Tomonari Furukawa.  Neural network constitutive modelling for non-linear characterization of anisotropic materials.  International journal for numerical methods in engineering, 85(8):939–957, 2011.
  * [75] Guanghui Liang and K Chandrashekhara.  Neural network based constitutive model for elastomeric foams.  Engineering structures, 30(7):2002–2011, 2008.
  * [76] V Papadopoulos, G Soimiris, DG Giovanis, and M Papadrakakis.  A neural network-based surrogate model for carbon nanotubes with geometric nonlinearities.  Computer Methods in Applied Mechanics and Engineering, 328:411–430, 2018.
  * [77] Daniel A White, William J Arrighi, Jun Kudo, and Seth E Watts.  Multiscale topology optimization using neural network surrogate models.  Computer Methods in Applied Mechanics and Engineering, 346:1118–1135, 2019.
  * [78] Liwei Wang, Yu-Chin Chan, Zhao Liu, Ping Zhu, and Wei Chen.  Data-driven metamaterial design with laplace-beltrami spectrum as “shape-dna”.  Structural and multidisciplinary optimization, 61:2613–2628, 2020.
  * [79] Liwei Wang, Yu-Chin Chan, Faez Ahmed, Zhao Liu, Ping Zhu, and Wei Chen.  Deep generative modeling for mechanistic-based learning and design of metamaterial systems.  Computer Methods in Applied Mechanics and Engineering, 372:113377, 2020.
  * [80] Zhaxylyk A Kudyshev, Alexander V Kildishev, Vladimir M Shalaev, and Alexandra Boltasseva.  Machine-learning-assisted metasurface design for high-efficiency thermal emitter optimization.  Applied Physics Reviews, 7(2), 2020.
  * [81] Doksoo Lee, Wei Chen, Liwei Wang, Yu-Chin Chan, and Wei Chen.  Data-driven design for metamaterials and multiscale systems: A review.  Advanced Materials, 36(8):2305254, 2024.
  * [82] Wei Ma, Feng Cheng, Yihao Xu, Qinlong Wen, and Yongmin Liu.  Probabilistic representation and inverse design of metamaterials based on a deep generative model with semi-supervised learning strategy.  Advanced Materials, 31(35):1901111, 2019.
  * [83] Hunter T Kollmann, Diab W Abueidda, Seid Koric, Erman Guleryuz, and Nahil A Sobh.  Deep learning for topology optimization of 2d metamaterials.  Materials & Design, 196:109098, 2020.
  * [84] Tianyu Zhao, Yiwen Li, Lei Zuo, and Kai Zhang.  Machine-learning optimized method for regional control of sound fields.  Extreme Mechanics Letters, 45:101297, 2021.
  * [85] Tianshuo Qiu, Xin Shi, Jiafu Wang, Yongfeng Li, Shaobo Qu, Qiang Cheng, Tiejun Cui, and Sai Sui.  Deep learning: a rapid and efficient route to automatic metasurface design.  Advanced Science, 6(12):1900128, 2019.
  * [86] Itzik Malkiel, Michael Mrejen, Achiya Nagler, Uri Arieli, Lior Wolf, and Haim Suchowski.  Plasmonic nanostructure design and characterization via deep learning.  Light: Science & Applications, 7(1):60, 2018.
  * [87] Rebekka V Woldseth, Niels Aage, J Andreas Bærentzen, and Ole Sigmund.  On the use of artificial neural networks in topology optimisation.  Structural and Multidisciplinary Optimization, 65(10):294, 2022.
  * [88] Aaditya Chandrasekhar and Krishnan Suresh.  Tounn: Topology optimization using neural networks.  Structural and Multidisciplinary Optimization, 63:1135–1149, 2021.
  * [89] Aaditya Chandrasekhar, Saketh Sridhara, and Krishnan Suresh.  Graded multiscale topology optimization using neural networks.  Advances in Engineering Software, 175:103359, 2023.
  * [90] Jorge Nocedal and Stephen J Wright.  Numerical optimization.  Springer, 1999.
  * [91] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [92] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [93] Aaditya Chandrasekhar, Saketh Sridhara, and Krishnan Suresh.  Auto: a framework for automatic differentiation in topology optimization.  Structural and Multidisciplinary Optimization, 64(6):4355–4365, 2021.
  * [94] Ivan Vasilev, Daniel Slater, Gianmario Spacagna, Peter Roelants, and Valentino Zocca.  Python Deep Learning: Exploring deep learning techniques and neural network architectures with Pytorch, Keras, and TensorFlow.  Packt Publishing Ltd, 2019.
  * [95] Erik Andreassen, Anders Clausen, Mattias Schevenels, Boyan S Lazarov, and Ole Sigmund.  Efficient topology optimization in matlab using 88 lines of code.  Structural and Multidisciplinary Optimization, 43:1–16, 2011.
  * [96] Jun Wu, Niels Aage, Rüdiger Westermann, and Ole Sigmund.  Infill optimization for additive manufacturing—approaching bone-like porous structures.  IEEE transactions on visualization and computer graphics, 24(2):1127–1140, 2017.
  * [97] Sarah N Hankins and Ray S Fertig III.  Methodology for optimizing composite design via biological pattern generation mechanisms.  Materials & Design, 197:109208, 2021.


# Hopfield-Enhanced Deep Neural Networks for Artifact-Resilient Brain State
Decoding

Arnau Marin-Llobet
Harvard University & FRCB-IDIBAPS
Allston, USA
<EMAIL_ADDRESS>
Arnau Manasanch
FRCB-IDIBAPS
Barcelona, Spain
<EMAIL_ADDRESS>
Maria V. Sanchez-Vives
FRCB-IDIBAPS & ICREA
Barcelona, Spain
<EMAIL_ADDRESS>

###### Abstract

The study of brain states, ranging from highly synchronous to asynchronous
neuronal patterns like the sleep-wake cycle, is fundamental for assessing the
brain’s spatiotemporal dynamics and their close connection to behavior.
However, the development of new techniques to accurately identify them still
remains a challenge, as these are often compromised by the presence of noise,
artifacts, and suboptimal recording quality. In this study, we propose a two-
stage computational framework combining Hopfield Networks for artifact data
preprocessing with Convolutional Neural Networks (CNNs) for classification of
brain states in rat neural recordings under different levels of anesthesia. To
evaluate the robustness of our framework, we deliberately introduced noise
artifacts into the neural recordings. We evaluated our hybrid Hopfield-CNN
pipeline by benchmarking it against two comparative models: a standalone CNN
handling the same noisy inputs, and another CNN trained and tested on
artifact-free data. Performance across various levels of data compression and
noise intensities showed that our framework can effectively mitigate
artifacts, allowing the model to reach parity with the clean-data CNN at lower
noise levels. Although this study mainly benefits small-scale experiments, the
findings highlight the necessity for advanced deep learning and Hopfield
Network models to improve scalability and robustness in diverse real-world
settings.

##  1 Introduction

Brain states are specific configurations of neural activity characterized by
distinct patterns of neuronal firing, spatiotemporal dynamics and behaviour
[1]. These configurations arise in various behavioral and cognitive conditions
and their understanding is fundamental to both neuroscience research and
clinical practice [2], since they underpin the brain’s ability to process
information and interact with the external world. This knowledge holds
significant implications for understanding both normal brain function and the
treatment of neurological disorders [3].

Natural brain states such as sleep versus wakefulness, or pharmacologically-
induced brain states such as under anesthetics provide functional models for
neuroscience research, as they enable the induction of specific spatiotemporal
cortical activity patterns [3]. In this context, deep anesthesia is
characterized by synchronous slow oscillations (SO) in the cortex [4]. These
SOs alternate between active ’Up states’ and quiescent ’Down states,’ offering
simpler models for neurological investigation. Here we will concentrate in
transitions involving SOs, wakefulness (AW), which shows desynchronized
activity, and microarousals (MA) that are brief, transient shifts towards
wakefulness that interrupt SO patterns. Understanding these global brain
states—AW, SO, and MA - and the transitions between them is still a key
challenge for neuroscience. A part of this challenge remains in the analysis
required to identify these states efficiently, particularly due to the need
for data compression strategies to handle large time-series datasets and
effective artifact removal techniques to negate experimental noise [5].

In this study, we employ Convolutional Neural Networks (CNNs) and pre-
processing via Hopfield Networks to improve the classification of artifact-
afflicted anesthesia brain states. Using low-dimensional binary images of
high-Signal-to-Noise Ratio (SNR) neural channels, we explore trade-offs
between accuracy, computational efficiency, and artifact distortion.

##  2 Related Work

Computational Methods for Brain Data. Computational advances have
revolutionized neural data analysis. While traditional statistical methods
remain important for basic research [6], machine learning techniques like CNNs
have gained prominence for tasks like anomaly detection and artifact-
contaminated signal decoding [5]. These methods also enable state
classification through interpretable deep learning [7]. On the clinical side,
breakthroughs have occurred in fields like robotics and brain-computers
interfaces [8], as well as cutting-edge neuroprosthetics for speech decoding
[9].

Hopfield Networks and NeuroAI Applications Hopfield Networks, originally
introduced by Hopfield in 1982 [10], have gained renewed attention for their
ability to act as associative memories in complex neural networks [11]. Modern
variants of these networks, often trained via backpropagation, have shown
superior memory storage properties and have been used as submodules in larger
AI networks [12]. Additionally, the interplay between associative memory
models and neuroscience has been a burgeoning area of research, offering bi-
directional insights that enrich both fields [13], [14]. Particularly relevant
to this study are applications of Hopfield Networks for denoising and
reconstructing neural data, showcasing the potential of these models to
enhance our understanding of complex neural systems [15][16].

##  3 Methodology

This study proposes a two-stage pipeline for classifying brain states under
varying anesthesia levels in rats. We selected the neural channel with the
highest SNR, transformed it into a lower-resolution binary image, and
introduced artifacts for evaluation of classification robustness. Artifact
correction is performed using a simple Hopfield Network, and final
classification is done using a CNN. See Fig.1.

Figure 1: Pipeline Overview. Top row illustrates the general workflow of our
methodology. Bottom row provides examples of various levels of artifact
addition within a single input image, as well as the key layers involved in
the model architecture for output generation.

###  3.1 Data Reconstruction and Artifact Addition

In our study, we analyze intracortical local field potentials (LFPs) from rats
under different anesthesia levels, sampled at 10 kHz with manual labeling for
each state. We segment the time-series data into 0.25-second intervals and
transform them into binary images through binarization. We then adjust the
image size by varying the compression ratio from $\frac{1}{25}$ to
$\frac{4}{25}$. After binarization, we introduce artifacts by altering pixel
values, guided by the err_percentage parameter. For a more comprehensive
explanation, see Appendix A.

###  3.2 Hopfield Network Model for Data Reconstruction

To address the inherent capacity constraints of Hopfield Networks, we employ a
multi-step approach. Firstly, a k-means clustering algorithm is used on the
binary images for each of the states of anesthesia to select representative
patterns. Each cluster’s centroid serves as a representative pattern, thus
allowing us to capture the underlying data structure while adhering to the
max_patterns limitation [10]. Next, individual Hopfield Networks are trained
for each neural state using Hebbian learning. Specifically, the Hebbian rule
is applied as $w_{ij}=\sum_{p}x_{i}^{(p)}x_{j}^{(p)}$, where $w_{ij}$ is the
weight between neurons $i$ and $j$, and $x^{(p)}$ is the $p^{th}$ pattern.
Note that for our implementation, self-connections $w_{ii}$ are set to zero.
These networks incorporate only artifact-added images, with the synaptic
weights stored in trained_weights. Finally, a template matching algorithm [17]
is deployed for data reconstruction, leveraging Hamming distance as a
similarity metric. The Hamming distance $D_{H}$ between two images $A$ and $B$
is computed as $D_{H}(A,B)=\sum_{i=1}^{n}|a_{i}-b_{i}|$, where $n$ is the
total number of pixels in the image, and $a_{i}$ and $b_{i}$ are individual
pixels in $A$ and $B$ respectively. This enables the selection of the most
similar clean image from the Hopfield Network for each artifact-corrupted
image, thus completing the reconstruction process.

###  3.3 Convolutional Neural Network (CNN)

The CNN takes artifact-corrected images as input, forming a supervised dataset
that is labeled according to the originating neural states. The CNN, tailored
for binary images, uses two convolutional and max pooling layers, followed by
two dense layers with ReLU activations, ending in a softmax classification
layer. Despite its smaller scale, the architecture balances resource
efficiency with task-specific robustness. This design choice is corroborated
by previous literature, for example as referenced in [18]. For more details on
model training and optimization, please refer to Appendix A.

###  3.4 Comparison and Validation of Results

Assessing the efficacy of our hybrid model for brain state classification is
complex due to an absence of standardized benchmarks. Traditional methods
often rely on manual curation and visual assessment, lacking quantitative
metrics for direct comparison. To tackle this, we compared the performance of
our combined Hopfield Network and CNN approach against two alternative
methodologies. To reduce variability, each model was run five times using
different random seeds.

The first alternative employs an identical CNN model, but omits both the
artifact introduction and the Hopfield Network reconstruction, feeding the CNN
directly with the binary images. The second alternative includes artifact
introduction but bypasses Hopfield Network-based reconstruction, providing the
CNN with artifact-corrupted images for classification. These comparisons allow
us to isolate the contributions of the artifact-handling elements in our
hybrid model and gauge their impact on classification accuracy.

##  4 Experimental Results

Figure 2: Experimental Results: The figure assesses classification accuracy
(test_acc) using three strategies: CNN-only without artifacts (shown in red),
Hopfield Network + CNN (in blue), and a baseline incorporating artifacts (in
green). Subfigure a) explores how test_acc varies with err_percentage,
aggregated over all image sizes. The accuracy of the CNN-only method without
artifacts stays consistent across various error percentages, as it repeatedly
processes identical inputs, akin to having different initialization seeds
despite a constant zero error rate. Subfigure b) illustrates the relationship
between test_acc and image size, keeping err_percentage range of 0.05 to 0.1
(highlighted in orange squares).

We employed Hopfield Networks for artifact removal in neural time series data,
aimed at classifying transitions in brain states from deep anesthesia to
wakefulness. This two-stage pipeline was benchmarked against two alternative
methods: one that utilizes a neural model without artifact addition, using
only clean signals; and a baseline strategy that retained artifacts but
omitted Hopfield Network pre-processing. Our evaluation focuses on the
efficacy of these various methods in brain state classification.

Figure 2.A indicates that at low err_percentage levels (0.05, 0.1), test_acc
is comparable between Hopfield Network reconstruction and CNN-only methods,
especially at an image size of 40 where original frequencies are better
preserved. At smaller image sizes (10, 20), less accurate frequency
representation leads to similar test_acc performance across all methods. For
err_percentage exceeding 0.2, the image distortions become too severe for
effective reconstruction, aligning the test_acc more closely with the
baseline, as observed in Figure 2.B.

##  5 Discussion and Future Work

We adopted a simple compression method for computational efficiency, at the
expense of model accuracy. Future research should investigate advanced
compression algorithms, such as autoencoders, to enhance state differentiation
and frequency representation. Fourier analysis might be a promising avenue
given its efficacy in capturing distinct frequencies characteristic of
different anesthesia states like AW and SO [19]. Moreover, refining the model
architecture could involve partitioning the MA states into asynchronous and
synchronous components, as suggested by recent literature [6].

This study offers useful insights for small-scale projects but warrants
further refinement for broader applications, particularly given challenges
like inter-subject variability. Our results suggest the efficacy of Hopfield
Networks for artifact reconstruction brain state classification under
anesthesia, yet these findings are preliminary. Future work should scrutinize
the robustness of this method across different distortion types and
compression levels. Despite constraints on the number of patterns we can
store, these limitations seem to have minimal impact on classification
accuracy in small-scale experiments. We posit that the high similarity within
short data segments of the same anesthesia state mitigates potential template-
matching errors. However, comprehensive evaluation is needed to quantify these
effects and uncover additional limitations and strengths.

##  6 Availability

All relevant code and demonstration dataset is publicly available at
https://github.com/arnaumarin/HDNN-ArtifactBrainState.

##  7 Acknowledgments and Disclosure of Funding

This work was made possible by European Union’s Horizon 2020 Framework
Programme for Research and Innovation under the Specific Grant Agreement No.
945539 (Human Brain Project SGA3) and by CORTICOMOD PID2020-112947RB-I00
funded by MCIN/ AEI /10.13039/501100011033.. We thank Leonardo Dalla Porta for
his support, Melody Torao for obtaining originally the data, and we thank and
acknowledge Cyprien Dcunha open source notebooks.

## References

  * [1] Harris, K. D., & Thiele, A. (2011). Cortical state and attention. Nature reviews neuroscience, 12(9), 509-523.
  * [2] Gervasoni, D., Lin, S. C., Ribeiro, S., Soares, E. S., Pantoja, J., & Nicolelis, M. A. (2004). Global forebrain dynamics predict rat behavioral states and their transitions. Journal of Neuroscience, 24(49), 11137-11147.
  * [3] Rosanova, M., Fecchio, M., Casarotto, S., Sarasso, S., Casali, A. G., Pigorini, A., … & Massimini, M. (2018). Sleep-like cortical OFF-periods disrupt causality and complexity in the brain of unresponsive wakefulness syndrome patients. Nature communications, 9(1), 4427.
  * [4] Sanchez-Vives, M. V., Massimini, M., & Mattia, M. (2017). Shaping the default activity pattern of the cortical network. Neuron, 94(5), 993-1001.
  * [5] Li, J., Struzik, Z., Zhang, L., & Cichocki, A. (2015). Feature learning from incomplete EEG with denoising autoencoder. Neurocomputing, 165, 23-31
  * [6] Tort-Colet, N., Capone, C., Sanchez-Vives, M. V., & Mattia, M. (2021). Attractor competition enriches cortical dynamics during awakening from anesthesia. Cell Reports, 35(12).
  * [7] Lee, M., Sanz, L. R., Barra, A., Wolff, A., Nieminen, J. O., Boly, M., … & Lee, S. W. (2022). Quantifying arousal and awareness in altered states of consciousness using interpretable deep learning. Nature communications, 13(1), 1064.
  * [8] Bonizzato, M., Pidpruzhnykova, G., DiGiovanna, J., Shkorbatova, P., Pavlova, N., Micera, S., & Courtine, G. (2018). Brain-controlled modulation of spinal circuits improves recovery from spinal cord injury. Nature communications, 9(1), 3015.
  * [9] Metzger, S. L., Littlejohn, K. T., Silva, A. B., Moses, D. A., Seaton, M. P., Wang, R., … & Chang, E. F. (2023). A high-performance neuroprosthesis for speech decoding and avatar control. Nature, 1-10.
  * [10] Hopfield, J. J. (1982). Neural networks and physical systems with emergent collective computational abilities. Proceedings of the national academy of sciences, 79(8), 2554-2558.
  * [11] Lucibello, C., & Mézard, M. (2023). The Exponential Capacity of Dense Associative Memories. arXiv preprint arXiv:2304.14964.
  * [12] Krotov, D., & amp; Hopfield, J. J. (2016). Dense associative memory for pattern recognition. Advances in Neural Information Processing Systems 29.
  * [13] Whittington, J. C., Warren, J., & Behrens, T. E. (2021). Relating transformers to models and neural representations of the hippocampal formation. arXiv preprint arXiv:2112.04035.
  * [14] Krotov, D., & Hopfield, J. (2020). Large associative memory problem in neurobiology and machine learning. arXiv preprint arXiv:2008.06996.
  * [15] Burns, T. F., & Fukai, T. (2023). Simplicial hopfield networks. arXiv preprint arXiv:2305.05179.
  * [16] Bricken, T., & Pehlevan, C. (2021). Attention approximates sparse distributed memory. Advances in Neural Information Processing Systems, 34, 15301-15315.
  * [17] Brunelli, R. (2009). Template matching techniques in computer vision: theory and practice. John Wiley & Sons.
  * [18] Tsinalis, O., Matthews, P. M., Guo, Y., & Zafeiriou, S. (2016). Automatic sleep stage scoring with single-channel EEG using convolutional neural networks.arXiv preprint arXiv:1610.01683.
  * [19] Torao-Angosto, M., Manasanch, A., Mattia, M., & Sanchez-Vives, M. V. (2021). Up and down states during slow oscillations in slow-wave sleep and different levels of anesthesia. Frontiers in systems neuroscience, 15, 609645.
  * [20] Islam, M. K., Rastegarnia, A., Nguyen, A. T., & Yang, Z. (2014). Artifact characterization and removal for in vivo neural recording. Journal of neuroscience methods, 226, 110-123.

##  8 Appendix A: Experimental Setup and Training Details

Experimental Dataset: Our demo-study focuses on approximately 5 minutes
extracted from a 6-hour intracortical local field potential (LFP) recording
from a rat under ketamine-medetomidine anesthesia. Each anesthesia state is
adequately represented. The dataset is sampled at a frequency of 10 kHz and
manually labeled by experts for ground truth annotations. The intracortical
LFP time-series is divided into 0.25-second segments, chosen to capture
stationary features across different anesthesia states. These segments are
transformed into binary images using a custom binarization algorithm pipeline.

In this pipeline, each 0.25-second segment is first normalized based on a
globally computed minimum and maximum value, ensuring consistent intensity
scaling across different states. Each normalized time-value is then mapped to
a pixel coordinate in the binary image; the x-coordinate corresponds to the
time index, and the y-coordinate to the normalized value. This transformation
effectively spatializes the time-series data into a 2D binary image, capturing
the waveform shapes inherent to the original series. To explore the impact of
compression, we manipulate the image_size between 10 and 40, corresponding to
compression ratios of 1/25 and 4/25, respectively. This image-based
representation serves as a practical and simple medium for neural data
analysis, allowing for varying degrees of compression.

After the binarization step, we introduce artificial artifacts by
algorithmically altering a predetermined subset of pixels. This subset is
determined by the err_percentage parameter, which serves to emulate common
neural recording noise and artifacts [20]. The err_percentage values can range
between 0 and 0.5, signifying that up to 50% of the total pixels may be
corrupted.

Ground truth annotations: Manual labeling for the establishment of ground
truth annotations was carried out by an expert who was instrumental in data
collection. This process involved defining three significant global states
stated before: AW, SO, and MA. Ground truth annotations were applied to the
LFP channel exhibiting the highest SNR. Additionally, an electromyography
(EMG) channel, positioned around the animal’s neck, was utilized to monitor
animal movements and served as a supplementary source of information for this
labeling process. This dataset originates from a series of experiments
designed to assess various brain states occurring within the rat cortical
brain under the influence of ketamine and medetomidine anesthesia. Each
experiment spanned approximately 5-6 hours, as shown in Figure 3. The entire
recording was categorized as follows.

Figure 3: Schematic of the different Ground Truth brain states. LFP of a full
transition is shown. It starts at Awake (orange) and after the anesthesia
induction it rapidly falls into the SO (green) regime. After a few hours it
goes to the MA state (blue), approaching wakefulness, as can be seen by the
EMG amplitude.

The AW state pertained to time periods showing asynchronous cortical activity
recorded before anesthesia administration. During this phase, the animal
exhibited clear signs of wakefulness, such as responsiveness and movement,
accompanied by a high EMG amplitude. SO labels were assigned to time segments
characterized by a distinct Up/Down pattern with a frequency below 2.5 Hz.
These segments immediately followed anesthesia induction and typically
persisted for 2-3 hours. The animal displayed no signs of wakefulness, and the
EMG exhibited a flat profile. MA labels were designated for time samples that
occurred once the stable SO state had concluded. These samples were expected
to exhibit a shift toward wakefulness, as detailed in [6]. They consisted of
brief periods of asynchronous activity interspersed with intervals of slow
oscillation characterized by a notably higher frequency (>2.5 Hz) compared to
the typical SO observed after anesthesia induction. During the MA state, the
animal remained unresponsive but began to display drowsiness-related symptoms,
which were validated through EMG and video recordings.

Any time samples showing noise or the presence of electrical or movement
artifacts were categorized as "artifact". Samples where definitive decisions
could not be reached were classified as "unknown".

CNN Model Training and Optimization:

The CNN takes artifact-corrected images as input, forming a supervised dataset
that is labeled according to the originating neural states. Engineered to
operate on binary images, the CNN architecture incorporates two 2D
convolutional layers with a kernel size of 3, each succeeded by max pooling
layers. The architecture concludes with two fully connected dense layers that
use ReLU activations, except for the softmax-activated final layer dedicated
to classification. It should be noted that although our architecture is
smaller compared to typical benchmarks in the field, it was designed to be
resource-efficient while still being sufficiently robust for the tasks at
hand.

The Convolutional Neural Network (CNN) model is subjected to a grid search to
optimize the number of epochs and the learning rate; the options considered
are [5, 10, 15], and [0.1, 0.001, 0.0001] respectfully. It is important to
note that we did not optimize the other hyperparameters because the primary
goal of this study is not to maximize and find the optimal classification
performance. Instead, we focus on assessing the Hopfield network’s ability to
reconstruct signals affected by artifacts.

Regarding dataset partitioning, 10% of the data is allocated for training, and
another 10% is set aside as a validation set to monitor the model and prevent
overfitting. The remaining 80% of the data serves as the test set. To evaluate
the effectiveness of the model, we use test accuracy, denoted as test_acc, and
a normalized confusion matrix, referred to as cm_normalized, as our key
performance indicators.