So far, the unknowns of our problem in Lagrangian approach have been the cell number density \(c_{t}\) , the cell deformation map \(y\)  and bound pili number density \(p_{0}\) .
A further challenge that arises in our problem is the presence of the Hessian of the spatial cell number density \(\nabla ^{2}_{x}c_{t}\)  and the Laplacian of the spatial cell number density \(\Delta _{x}c_{t}\)  in the bound pili number density evolution equation (REF ) and the pili density formation Piola stress equation (REF ).
Thus, second derivatives of the spatial cell number density \(c_t\)  are required to be calculated.
In doing so, different strategies have been employed in the literature, among which the well-established ones are \(\mathcal {C}^{1}\)  continuous elements [1]}, [2]}, isogeometric analysis [3]}, [4]}, micromorphic continuum approach and the gradient enhanced framework [5]}, [6]}, [7]}.
In order to stay within the realm of the classical finite element method (FEM) associated with a \(\mathcal {C}^0\)  continuous interpolation approach, we adopt the gradient enhanced framework.
In doing so, we introduce an additional independent spatial vector field \(g\)  to represent the spatial gradient of the cell number density \(\nabla _{x} c_t\)  and weakly enforce \(g\overset{!}{=}\nabla _{x} c_t\) .
Note, similar to the spatial cell number density, its gradient is also parametrized as a function of the Lagrangian coordinates \(g=g(X,t)\) .
The first step is to transform the Laplacian and Hessian to substitute the terms including spatial derivatives of \(c_{t}\)  with \(g\)  as
\(\begin{aligned}&\Delta _{x} c_t = \nabla _{x} \cdot (\nabla _{x} c_t) = \nabla _{x} \cdot g \,,\\&\nabla ^{2}_{x} c_t = \nabla _{x} (\nabla _{x} c_t) = \nabla _{x}^{\text{sym}} g = {2} \left[\nabla _{x} g + \nabla _{x}^{T} g \right]\,,\end{aligned}\) 
