paper_vem_section.tex

% =============================================================================
% Section 5: Virtual Element Method for Biofilm Mechanics
% Draft for inclusion in the main paper (Nishioka 2026, IKM Hannover)
% =============================================================================

\section{Future Directions: Virtual Element Method for Biofilm Mechanics}
\label{sec:vem}

The finite element analyses presented in Sections~3--4 rely on structured meshes
(hexahedral Q4 in~2D, tetrahedral C3D4 in~3D) generated from confocal images
through a 5-step pipeline: confocal $\to$ voxel $\to$ marching cubes $\to$
tetrahedral mesh $\to$ Abaqus input.
In this section we present a Virtual Element Method (VEM) framework that
reduces this to a \emph{2-step pipeline} (confocal $\to$ Voronoi $\to$ VEM)
while naturally accommodating the features most relevant to biofilm mechanics:
arbitrary element shapes, per-colony material assignment, viscoelastic
constitutive response, and adaptive refinement.

% ---------------------------------------------------------------------------
\subsection{VEM formulation for linear elasticity}
\label{sec:vem:elasticity}

The VEM~\cite{BeiraodaVeiga2013,BeiraodaVeiga2014} solves elliptic PDEs on
meshes of \emph{arbitrary polygonal or polyhedral elements} without
constructing explicit shape functions.  For lowest-order ($k{=}1$) 2D
elasticity, the local virtual element space on a polygon~$E$ with $n_v$
vertices is
\begin{equation}
  V_h^E = \bigl\{\, \mathbf{v}_h \in [H^1(E)]^2 :
          \mathbf{v}_h|_{\partial E} \text{ piecewise linear},\,
          \Delta \mathbf{v}_h = \mathbf{0} \text{ in } E \,\bigr\},
\end{equation}
with degrees of freedom $\{\mathbf{u}_i\}_{i=1}^{n_v}$ at the vertices,
giving $n_{\mathrm{dof}}^E = 2n_v$ per element.

The key construction is the \emph{elliptic projection}
$\Pi^\nabla : V_h^E \to [\mathcal{P}_1(E)]^2$
that maps virtual functions onto the 6-dimensional space of linear
polynomials (3 rigid-body modes $+$ 3 constant strain modes).
The element stiffness matrix decomposes as
\begin{equation}
  \mathbf{K}^E = \mathbf{K}^E_\pi + \mathbf{K}^E_\mathrm{stab},
  \qquad
  \mathbf{K}^E_\pi = (\Pi^\nabla)^\mathsf{T}\, \tilde{\mathbf{C}}\, \Pi^\nabla,
  \label{eq:vem_stiffness}
\end{equation}
where $\tilde{\mathbf{C}}$ incorporates the plane-stress constitutive
matrix~$\mathbf{C}$ and element area~$|E|$, and
$\mathbf{K}^E_\mathrm{stab} = \alpha_s\, \mathrm{tr}(\mathbf{C})\,|E|\,
(\mathbf{I} - \Pi^\nabla\,\mathbf{D})^\mathsf{T}
(\mathbf{I} - \Pi^\nabla\,\mathbf{D})$
stabilises the kernel of the projection with $\alpha_s = 0.5$.

For 3D polyhedra (Section~\ref{sec:vem:3d}), the polynomial space extends
to 12~modes (3 translations $+$ 3 rotations $+$ 6 constant strain),
with face-based integration replacing edge-based integration in the
$\mathbf{B}$ and $\mathbf{D}$ matrices~\cite{Wriggers2024}.

\paragraph{Verification.}
A manufactured-solution convergence study on Voronoi meshes yields
$L^2$ rate~$\approx 2.14$ and $H^1$ rate~$\approx 1.29$, exceeding the
theoretical minimum ($L^2 \geq 2$, $H^1 \geq 1$) and outperforming
linear triangular FEM ($L^2 = 1.88$, $H^1 = 0.99$) on the same problem
(Table~\ref{tab:vem_convergence}).

\begin{table}[h]
\centering
\caption{$h$-convergence rates: VEM vs FEM on Cook's membrane problem.}
\label{tab:vem_convergence}
\begin{tabular}{lcc}
\toprule
Method & $L^2$ rate & $H^1$ rate \\
\midrule
VEM (Voronoi) & 2.14 & 1.29 \\
VEM (quadrilateral) & 2.03 & 1.99 \\
FEM (triangle) & 1.88 & 0.99 \\
\bottomrule
\end{tabular}
\end{table}

% ---------------------------------------------------------------------------
\subsection{Viscoelastic VEM (VE-VEM) with Simo~1987 integrator}
\label{sec:vem:viscoelastic}

Biofilm is a viscoelastic material~\cite{Peterson2015,Flemming2010} whose
relaxation time depends on species composition.  We model this using a
Standard Linear Solid (SLS / Zener) with DI-dependent parameters:
\begin{equation}
  E_\infty(\mathrm{DI}) = E_{\min} + (E_{\max}-E_{\min})(1-\mathrm{DI})^n,
  \quad
  E_0 = 2\,E_\infty,
  \quad
  \tau(\mathrm{DI}) = \tau_{\min} + (\tau_{\max}-\tau_{\min})(1-\mathrm{DI})^n,
  \label{eq:sls_params}
\end{equation}
with $E_{\max}=1000$\,Pa, $E_{\min}=10$\,Pa, $\tau_{\max}=60$\,s,
$\tau_{\min}=2$\,s, $n=2$.  The viscous branch modulus is
$E_1 = E_0 - E_\infty$ and viscosity $\eta = E_1\,\tau$.

The constitutive relation in Voigt form reads
\begin{equation}
  \boldsymbol{\sigma}(t)
  = \mathbf{C}_\infty\,\boldsymbol{\varepsilon}(t) + \mathbf{h}(t),
  \label{eq:sls_stress}
\end{equation}
where $\mathbf{h}$ is the internal (viscous) stress variable governed by the
evolution ODE
$\dot{\mathbf{h}} = \mathbf{C}_1\,\dot{\boldsymbol{\varepsilon}}
- \mathbf{h}/\tau$.
Following Simo~\cite{Simo1987}, we integrate this exactly over a time step
$[t_n, t_{n+1}]$ with $\Delta t = t_{n+1} - t_n$:
\begin{equation}
  \mathbf{h}_{n+1}
  = e^{-\Delta t/\tau}\,\mathbf{h}_n
    + \gamma\,\mathbf{C}_1
      \bigl(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n\bigr),
  \quad
  \gamma = \frac{\tau}{\Delta t}\,\bigl(1 - e^{-\Delta t/\tau}\bigr).
  \label{eq:simo_update}
\end{equation}
The algorithmic tangent modulus is
$\mathbf{C}_\mathrm{alg} = \mathbf{C}_\infty + \gamma\,\mathbf{C}_1$,
leading to the time-step stiffness matrix
\begin{equation}
  \mathbf{K}_\mathrm{alg}^E
  = \bigl(\boldsymbol{\Pi}_\varepsilon^E\bigr)^\mathsf{T}\,
    \mathbf{C}_\mathrm{alg}^E\, |E|\,
    \boldsymbol{\Pi}_\varepsilon^E
  + \mathbf{K}_\mathrm{stab}^E,
  \label{eq:ve_vem_K}
\end{equation}
where $\boldsymbol{\Pi}_\varepsilon^E$ is the strain projector
($3 \times 2n_v$ in 2D, $6 \times 3n_v$ in 3D).

\paragraph{Validation.}
Under laterally confined uniaxial step strain
($u_x = 0$ for all nodes, $u_y = \varepsilon_0$ on top),
the analytical stress response is
\begin{equation}
  \sigma_{yy}(t)
  = \frac{E_\infty + E_1\,e^{-t/\tau}}{1-\nu^2}\,\varepsilon_0.
  \label{eq:confined_analytical}
\end{equation}
Our implementation reproduces this to \emph{machine precision}
in both 2D (max relative error $1.3 \times 10^{-15}$, 64 Voronoi cells)
and 3D (max relative error $4.9 \times 10^{-16}$, 27 hexahedral cells),
confirming the correctness of the Simo integrator and VEM projector
coupling.

% ---------------------------------------------------------------------------
\subsection{Growth-coupled VE-VEM}
\label{sec:vem:growth}

The staggered coupling of Section~2.3 extends naturally to the viscoelastic
setting.  At each growth step, the Hamilton 5-species ODE yields
time-dependent species fractions $\phi_i(t)$, from which we compute
\begin{equation}
  \mathrm{DI}(t) = -\sum_{i=1}^{5} \phi_i(t)\,\ln\phi_i(t) \;/\; \ln 5,
\end{equation}
and update the per-element SLS parameters via Eq.~\eqref{eq:sls_params}.
The VE-VEM then computes the viscoelastic response for the current time
interval $[t_n, t_{n+1}]$.

The growth-coupled solver reveals a key asymmetry:
\begin{itemize}
  \item \textbf{Commensal}: DI decreases ($0.38 \to 0.23$) as the ecosystem
    diversifies $\to$ $E_\infty$ increases ($481 \to 913$\,Pa),
    $\tau$ increases ($26 \to 41$\,s) $\to$ biofilm \emph{stiffens and relaxes
    more slowly}.
  \item \textbf{Dysbiotic}: DI increases ($0.40 \to 0.49$) as \emph{S.~oralis}
    dominates $\to$ $E_\infty$ decreases ($481 \to 324$\,Pa),
    $\tau$ decreases ($26 \to 23$\,s) $\to$ biofilm \emph{softens and relaxes
    faster}.
\end{itemize}
This divergent mechanical evolution under identical initial conditions
demonstrates that species composition dynamics directly drive the
time-dependent mechanical response, validating the DI-based constitutive
model in the viscoelastic regime.

% ---------------------------------------------------------------------------
\subsection{3D polyhedral VE-VEM}
\label{sec:vem:3d}

The 2D framework extends to 3D polyhedral elements following
Wriggers et al.~\cite{Wriggers2019,Wriggers2024}.  On a polyhedron~$E$ with
$n_v$ vertices and $n_f$ faces:
\begin{itemize}
  \item The polynomial basis has 12 modes
    ($\mathbf{p}_1, \ldots, \mathbf{p}_{12}$):
    3~translations, 3~rotations, 6~constant strain modes.
  \item The $\mathbf{B}$ matrix ($12 \times 3n_v$) is assembled via
    face integrals using the midpoint rule on triangulated faces.
  \item The strain projector $\boldsymbol{\Pi}_\varepsilon$ is
    $6 \times 3n_v$, extracted from rows 7--12 of the full projector.
  \item The algorithmic tangent uses $6 \times 6$ Voigt matrices
    $\mathbf{C}_\infty^{3\mathrm{D}}, \mathbf{C}_1^{3\mathrm{D}}$
    for the isotropic SLS model.
\end{itemize}
The 3D implementation achieves the same machine-precision validation
($4.9 \times 10^{-16}$) under confined compression, confirming
dimension-independent correctness of the Simo integrator.

% ---------------------------------------------------------------------------
\subsection{Additional VEM modules}
\label{sec:vem:additional}

Beyond the VE-VEM core, the framework includes:
\begin{itemize}
  \item \textbf{Neo-Hookean VEM} (A1): Large-deformation hyperelasticity
    with Newton-Raphson solver, demonstrating 43\% displacement difference
    versus linear theory at 1\% nominal strain.  For biofilm strains
    $\varepsilon > 5$--$10\%$ (typical under gingival crevicular fluid
    pressure), nonlinear effects are significant.
  \item \textbf{Phase-field fracture} (B1): Aldakheel et al.~\cite{Aldakheel2018}
    framework with DI-dependent fracture toughness
    $G_c(\mathrm{DI})$.  Dysbiotic regions ($G_c = 0.01\,\text{J/m}^2$)
    crack catastrophically while commensal periphery
    ($G_c = 0.5\,\text{J/m}^2$) remains intact.
  \item \textbf{Cohesive zone model} (B2): Bilinear traction-separation law
    for the tooth--biofilm interface, with DI-dependent cohesive strength
    $\sigma_{\max}$.
  \item \textbf{Adaptive refinement} (B3): A posteriori error estimator
    with Dörfler marking and automatic $h$-refinement at crack tips,
    demonstrated with 40 $\to$ 121~cells over 3~refinement levels.
  \item \textbf{P\textsubscript{2} VEM} (A2): Second-order elements with
    vertex $+$ edge-midpoint DOFs, achieving 15--45\% improvement in stress
    accuracy over P\textsubscript{1}.
  \item \textbf{Confocal $\to$ VEM pipeline}: Direct conversion of 5-channel
    fluorescence images (matching Heine~2025 species distributions) to
    Voronoi meshes with per-colony DI assignment.
\end{itemize}

The full framework comprises 12~solver modules validated by
120$+$ automated tests (Table~\ref{tab:vem_tests}).

\begin{table}[h]
\centering
\caption{Test coverage of the VEM framework.}
\label{tab:vem_tests}
\begin{tabular}{lrl}
\toprule
Module & Tests & Key validation \\
\midrule
Poisson (scalar) & 6 & Patch test, $h$-convergence \\
2D elasticity & 9 & Patch test ($10^{-18}$), manufactured solution \\
3D elasticity & 15 & Patch test ($10^{-19}$), Voronoi polyhedra \\
VE-VEM (2D) & 25 & Machine precision relaxation, DI ordering \\
P\textsubscript{2} VEM & 23 & Projector identity, analytical strain energy \\
Space-time VEM & 4 & Anisotropic Voronoi, SLS time stepping \\
Growth-coupled & 12 & Replicator ODE, cell division \\
Error estimator & 12 & ZZ-type, Dörfler marking \\
\midrule
\textbf{Total} & \textbf{106+} & \\
\bottomrule
\end{tabular}
\end{table}

% ---------------------------------------------------------------------------
\subsection{Significance for biofilm mechanics}
\label{sec:vem:significance}

The VEM framework addresses three fundamental limitations of conventional
FEM for biofilm analysis:

\begin{enumerate}
  \item \textbf{Mesh--image compatibility}: Confocal colony segmentation
    directly yields Voronoi cells that serve as VEM elements, eliminating
    the mesh generation bottleneck (2-step vs.\ 5-step pipeline).
  \item \textbf{Per-colony constitutive assignment}: Each Voronoi cell
    maps 1:1 to a micro-colony with its own DI, enabling spatially resolved
    DI $\to$ $E_\infty(\mathrm{DI})$, $\tau(\mathrm{DI})$, $G_c(\mathrm{DI})$
    assignment without sub-element averaging.
  \item \textbf{Topology changes}: Biofilm growth, detachment, and crack
    propagation alter the domain topology.  VEM's tolerance of arbitrary
    element shapes means that cell division, adaptive refinement, and
    crack insertion require no global remeshing.
\end{enumerate}

To the best of our knowledge, this is the first application of the
Virtual Element Method to biofilm mechanics, combining the VEM solid
mechanics expertise of the IKM Hannover group~\cite{Wriggers2019,
Wriggers2024,Aldakheel2018} with TMCMC-calibrated species dynamics.