A Finite Element Method for Modeling Diffusion and Drug Release from Nanocellulose/Nanoporous Silicon Composites

Abstract

Background and Objective: A previous study investigated the in vitro release of methylene blue (MB), a widely used cationic dye in biomedical applications, from nanocellulose/nanoporous silicon (NC/nPSi) composites under conditions simulating body fluids. The results showed that MB release rates varied significantly with the nPSi concentration in the composite, highlighting its potential for controlled drug delivery. To further analyze the relationship between diffusion dynamics and the MB concentration, this study developed a finite element (FE) method to solve Fick’s equations governing the drug delivery system. Methods: Release profiles of MB from NC/nPSi composites with varying nPSi concentrations (0%, 0.1%, 0.5%, and 1.0%) were experimentally measured in triplicate using phosphate-buffered saline (PBS) at 37 °C, pH 7.4, and 100 rpm. Mathematical models incorporating linear and quadratic dependencies of the diffusion coefficient on the MB concentration were developed and tested using the FE method. Model parameters were refined by minimizing the error between simulated and experimental MB release profiles. Results: The proposed FE method closely matched experimental data, validating its accuracy and robustness in simulating the diffusion and release processes. Conclusions: This study emphasizes the significant impact of the nPSi concentration on enhancing release control and highlights the importance of material composition in designing drug delivery systems. The findings suggest that the FE method can be effectively applied to model other complex systems, paving the way for advancements in precision drug delivery and broader biomedical applications.

1. Introduction

Nanoporous silicon (nPSi) stands out as an excellent biomaterial for drug delivery applications due to its high surface area, biocompatibility, biodegradability, and bioresorbability [1,2,3]. Typically, drugs are either loaded into the porous matrix or immobilized on the surface following appropriate surface derivatization. When combined with biopolymers, it acts as a substrate for composite materials, introducing advantageous chemical and physical properties not present in individual components. These benefits include improved control over drug release kinetics and enhanced stability in aqueous solutions [4,5]. As a result, nPSi has been paired with various organic matrices to create composites as advanced drug delivery systems, such as $\beta$-cyclodextrin polymers [4,6,7], oxidized hyaluronic acid hydrogels [8], and poly(L-lactide) acid [9], among others.

In this context, nanocellulose (NC) has recently emerged as one of the most promising “green” materials for obtaining drug delivery carriers as composites. It offers adaptable surface chemistry, a high surface area, biocompatibility, and biodegradability [10,11]. Recently, K. Garrido-Miranda et al. [12] synthesized NC/nPSi composites for the controlled release of methylene blue (MB), a cationic thiazine dye widely used in biomedicine for various purposes. These include the treatment of methemoglobinemia [13], its use as a marker and indicator in various surgical techniques [14], and its use as an analgesic in different treatments [15]. Additionally, it has been applied as an antibacterial, antiviral agent, and against cancer cells [16].

On the other hand, the finite element (FE) method is a numerical technique commonly used to approximate the solution of ordinary and partial differential equations (see, e.g., [17,18,19] and the references therein). It involves writing the equation in its variational (weak) form and approximating the solution, originally defined in an infinite-dimensional space, by restricting it to a finite-dimensional subspace. This approach has been successfully applied to simulate various physical phenomena, including the fluid–structure interaction [20], electromagnetism [21], and heat transfer [22].

This work presents an FE method to simulate the diffusion and controlled release of MB from NC/nPSi composites. The simulations are based on Fick’s second law [23], a widely used model for release kinetics (e.g., [24,25,26,27,28]). In this study, a concentration-dependent diffusion coefficient is incorporated to more accurately represent the behavior of highly polymerized NC. Furthermore, the influence of temperature on diffusivity in biological environments is considered a critical factor. This temperature dependency is modeled using the Stokes–Einstein relationship [29], offering a deeper insight into the dynamics of the diffusion process:

$$\begin{array}{c}\hfill D={\displaystyle \frac{k_{B}T}{a\eta r}},\end{array}$$

where $k_B$ is Boltzmann’s constant, T represents the temperature, a denotes the viscosity, $\eta$ is a constant that accounts for the boundary conditions between the diffusing molecule and the solvent, and r is the radius of the molecule. Consequently, diffusion coefficients increase with temperature, as observed in drug delivery experiments [30,31].

Nevertheless, as this model is further refined by incorporating dependency relations from [32], which align with the experimental release profiles reported by K. Garrido-Miranda et al. [12], the influence of temperature was not investigated and remains limited to a concentration-dependent diffusion coefficient.

Fick’s second law is first discretized using the FE method, transforming the governing equation into a discrete system. Since the diffusion coefficient depends on concentration, the non-linearities are addressed using the Picard iteration. This method linearizes the system at each time step, enabling its resolution using Octave [33]. The proof of unique solvability of the system is provided under the assumption of a bounded diffusion coefficient. Based on this assumption, examples are presented to offer valuable insights into designing more efficient drug delivery systems for biomedical applications.

The outline of this article is as follows. Section 2 presents the model problem and its FE discretization. A case study comparing the numerical solution with experimental data is presented in Section 3. Numerical details concerning future directions are discussed in Section 4, while conclusions are drawn in Section 5.

2. Materials and Methods

2.1. Model Problem

This work is about the one-dimensional diffusion and controlled release of MB through an NC/nPSi composite of thickness $2L$ (see Figure 1). It is assumed that L is sufficiently small, such that the amount of MB passing through the edges of the composite is negligible.

In this context, the drug delivery system is modeled within the framework of Fick’s second law [23], a widely used approach in diffusion processes, with a diffusion coefficient D that depends on the concentration $\varphi$. Let $I=(-L,L)$ be an interval along the x-axis, and let $T>0$ be a given finite time. The governing equation is

$${\displaystyle \frac{\partial \varphi }{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(D\left(\varphi \right){\displaystyle \frac{\partial \varphi }{\partial x}}\right)\mathrm{in}I\times (0,T).$$

(1)

At the surfaces $x=\pm L$, the concentration is prescribed as $\varphi ={\varphi }_{\infty }$ for all $t\in (0,T)$.

When the concentration is initially uniform, with $\varphi (x,0)={\varphi }_0$ for all $x\in I$, and both D and ${\varphi }_{\infty }$ are constants, the partial differential Equation (1) can be solved using the method of separation of variables or the Laplace transform, as detailed in [34]. This gives

$${\displaystyle \frac{\varphi (x,t)-{\varphi }_0}{{\varphi }_{\infty }-{\varphi }_0}}=1-{\displaystyle \frac{4}{\pi }}\sum _{n=0}^{\infty }{\displaystyle \frac{{(-1)}^n}{(2n+1)}}cos\left[{\displaystyle \frac{(2n+1)\pi x}{2L}}\right]exp\left[-{\displaystyle \frac{D{(2n+1)}^2{\pi }^{2}t}{4L^2}}\right].$$

(2)

Moreover, if $M_t$ denotes the amount of MB released at time t and $M_{\infty }$ the corresponding quantity after infinite time, the release profile is given by

$${\displaystyle \frac{M_t}{M_{\infty }}}=1-{\displaystyle \frac{8}{{\pi }^2}}\sum _{n=0}^{\infty }{\displaystyle \frac{1}{{(2n+1)}^2}}exp\left[-{\displaystyle \frac{D{(2n+1)}^2{\pi }^{2}t}{4L^2}}\right].$$

(3)

Although assuming a constant D is common in diffusion modeling, this assumption is unrealistic for highly polymerized substances [34], such as NC and NC-based composites. Given the relevance of the NC in this context, Fick’s second law in (1) with a variable D will be addressed in the following sections, focusing on the numerical solution. To facilitate understanding, we first analyze the case with a constant D and then extend the analysis to accommodate $D=D\left(\varphi \right)$.

2.2. Constant Diffusion Coefficient

2.2.1. Preliminaries

Let us review some basic concepts of functional analysis which are useful in dealing with partial differential equations. We first define the Sobolev space $H^1\left(I\right)$ as

$$H^1\left(I\right):=\{\psi \in L^2\left(I\right):{\psi }^{\prime }\in L^2\left(I\right)\},$$

which is a Hilbert space equipped with the inner product

$${(\varphi ,\psi )}_{1,I}:={\int }_I\varphi \psi \mathrm{d}x+{\int }_I{\varphi }^{\prime }{\psi }^{\prime }\mathrm{d}x\forall \varphi ,\psi \in H^1\left(I\right).$$

For further details, we refer the reader to [35]. The norm induced by ${(·,·)}_{1,I}$ is given by

$${\parallel \psi \parallel }_{1,I}:={\left({\parallel \psi \parallel }_{0,I}^2+{\left|\psi \right|}_{1,I}^2\right)}^{1/2}\forall \psi \in H^1\left(I\right),$$

where ${\left|\psi \right|}_{1,I}:={\parallel {\psi }^{\prime }\parallel }_{0,I}$ is a semi-norm on $H^1\left(I\right)$. It is well known that the quantity ${\left|\psi \right|}_{1,I}$ is a norm equivalent to ${\parallel \psi \parallel }_{1,I}$ on the following subspace of $H^1\left(I\right)$:

$$H_0^1\left(I\right):=\left\{\psi \in H^1\left(I\right):\psi =0\mathrm{at}x=\pm L\right\},$$

(4)

due to Poincaré’s inequality. This result is stated below, and its proof can be found in ([35], Proposition 8.13).

Lemma 1

(Poincaré’s inequality). Let I be a bounded interval. Then, there exists a constant $C_P>0$, depending only on I, such that

$${\parallel \psi \parallel }_{1,I}\le C_P{\left|\psi \right|}_{1,I}\forall \psi \in H_0^1\left(I\right).$$

(5)

Finally, we note that an alternative way of interpreting $\varphi$ in (1) is to treat it as a function of time, taking values in a Sobolev space, e.g., V, where the elements of V are functions that depend only on the spatial variable:

$$\varphi :t\in (0,T)\mapsto \varphi \left(t\right)\equiv \varphi (·,t)\in V.$$

This notation will be used throughout this work. Furthermore, time derivatives will be denoted by $d_t\left(·\right)$.

2.2.2. Weak Formulation

To introduce the weak formulation of Fick’s second law, we multiply (1) by a test function $\psi \in H_0^1\left(I\right)$, assume a constant D, and perform integration over I, yielding

$${\int }_{-L}^L\psi d_t\varphi \left(t\right)\mathrm{d}x+D{\int }_{-L}^L{\varphi }^{\prime }\left(t\right){\psi }^{\prime }\mathrm{d}x-{D{\varphi }^{\prime }\left(t\right)\psi \mid }_{x=-L}^{x=L}=0.$$

Then, using the condition that $\psi =0$ at $x=\pm L$ (cf. (4)), we obtain

$${\int }_{-L}^L\psi d_t\varphi \left(t\right)\mathrm{d}x+D{\int }_{-L}^L{\varphi }^{\prime }\left(t\right){\psi }^{\prime }\mathrm{d}x=0.$$

(6)

By virtue of identity (6) and under the conditions used to determine the exact $\varphi$ in (2), the weak formulation of Equation (1) with a constant D reads the following: For almost every $t\in (0,T)$, find $\varphi \left(t\right)\in H^1\left(I\right)$, such that $\varphi ={\varphi }_{\infty }$ on $\{-L,L\}\times (0,T)$, $\varphi ={\varphi }_0$ on $I\times \left\{0\right\}$, and 

$${\int }_{-L}^L\psi d_t\varphi \left(t\right)\mathrm{d}x+D{\int }_{-L}^L{\varphi }^{\prime }\left(t\right){\psi }^{\prime }\mathrm{d}x=0\forall \psi \in H_0^1\left(I\right).$$

(7)

At this point, we emphasize that the FE method, which will be detailed later, does not directly allow us to deduce the unique solvability for the discretization of (7) when ${\varphi }_{\infty }\ne 0$. This is because, in this case, the discrete versions of $\varphi \left(t\right)$ and $\psi$ reside in different spaces. To address this issue, we recall that the orthogonal complement of $H_0^1\left(I\right)$ in $H^1\left(I\right)$ is defined by

$$H_0^1{\left(I\right)}^{\perp }:=\{\lambda \in H^1\left(I\right):{(\lambda ,\psi )}_{1,I}=0\forall \psi \in H_0^1\left(I\right)\}.$$

(8)

It follows that $\lambda \in H_0^1{\left(I\right)}^{\perp }$ if and only if $\lambda$ is the weak solution to the equation

$$-{\lambda }^{\prime \prime }+\lambda =0.$$

Accordingly, we decompose $H^1\left(I\right)$ as the direct sum $H_0^1\left(I\right)\oplus W$, where W is the subspace spanned by $\{e^x,e^{-x}\}$. Next, we introduce the auxiliary unknown

$$\vartheta \left(t\right):=\varphi \left(t\right)-\lambda ,$$

(9)

with $\lambda \in H_0^1{\left(I\right)}^{\perp }$ being defined by

$$\lambda \left(x\right):=({\displaystyle \frac{{\varphi }_{\infty }}{e^L+e^{-L}}})(e^x+e^{-x}).$$

(10)

It is easy to check that $\lambda (-L)=\lambda \left(L\right)={\varphi }_{\infty }$, from which we deduce that (7) is equivalent to the following problem: for almost every $t\in (0,T)$, find $\vartheta \left(t\right)\in H_0^1\left(I\right)$ such that $\vartheta ={\varphi }_0$ on $I\times \left\{0\right\}$ and

$${\int }_{-L}^L\psi d_t\vartheta \left(t\right)\mathrm{d}x+D{\int }_{-L}^L{\vartheta }^{\prime }\left(t\right){\psi }^{\prime }\mathrm{d}x=-D{\int }_{-L}^L{\lambda }^{\prime }{\psi }^{\prime }\mathrm{d}x\forall \psi \in H_0^1\left(I\right).$$

(11)

We can therefore use $\vartheta$ to recover the solution to (7). In particular, $\vartheta =\varphi$ when ${\varphi }_{\infty }=0$.

2.3. Discretization of the Model with a Constant Diffusion Coefficient

This section focuses on approximating the solution to (11) using a fully discrete scheme that combines an FE method in space with an implicit Euler method in time. For simplicity, homogeneous boundary conditions are initially assumed. The non-homogeneous case will be addressed in Section 2.3.3.

2.3.1. Fully Discrete Scheme

Let $V_h$ denote an arbitrary finite-dimensional subspace of $H_0^1\left(I\right)$. We begin by considering the problem given by (11) with ${\varphi }_{\infty }=0$. We discretize this problem in space using the following FE scheme: For each $t\in [0,T]$, find ${\vartheta }_h\left(t\right)\in V_h$ such that ${\vartheta }_h^0={\varphi }_h^0$ and

$${\int }_{-L}^L{\psi }_h{d}_t{\vartheta }_h\left(t\right)\mathrm{d}x+A(t,{\vartheta }_h,{\psi }_h)=0\forall {\psi }_h\in V_h,$$

(12)

where the bilinear form $A:t\mapsto H_0^1\left(I\right)\times H_0^1\left(I\right)$ is defined by

$$A(t,{\vartheta }_h,{\psi }_h):=D{\int }_{-L}^L{\vartheta }_h^{\prime }\left(t\right){\psi }_h^{\prime }\mathrm{d}x,$$

(13)

and ${\varphi }_h^0$ is the $L^2$-projection of ${\varphi }_0$ into $V_h$, such that ${\varphi }_h^0\in V_h$ satisfies, for all ${\psi }_h\in V_h$,

$${\int }_{-L}^L{\varphi }_h^0{\psi }_h\mathrm{d}x={\int }_{-L}^L{\varphi }_0{\psi }_h\mathrm{d}x.$$

(14)

To discretize in time, we partition the interval $[0,T]$ as

$$0=t^0<t^1<\cdots <t^{N+1}=T,$$

with time step denoted by $\mathsf{\Delta }{t}^n:=t^{n+1}-t^n$ for $n\in \{0,1,\dots ,N\}$. Furthermore, we denote a function $\zeta \left(t\right)$ at time level $t=t^n$ by ${\zeta }^n$.

For the implicit Euler method in time, we use

$${\left(d_t{\vartheta }_h\right)}^{n+1}\approx {\displaystyle \frac{{\vartheta }_h^{n+1}-{\vartheta }_h^n}{\mathsf{\Delta }{t}^n}}.$$

Inserting this expression into (12) at time $t^{n+1}$ yields the following fully discrete formulation of (11) with homogeneous boundary conditions: For each $n\in \{0,1,\dots ,N\}$, find ${\vartheta }_h^{n+1}\in V_h$, such that

$$B({\vartheta }_h^{n+1},{\psi }_h)+\mathsf{\Delta }{t}^{n}A({\vartheta }_h^{n+1},{\psi }_h)=B({\vartheta }_h^n,{\psi }_h)\forall {\psi }_h\in V_h,$$

(15)

where, for easy of notation, we write $A({\vartheta }_h^{n+1},{\psi }_h)$ instead of $A(t^{n+1},{\vartheta }_h^{n+1},{\psi }_h)$ and

$$B({\vartheta }_h,{\psi }_h):={\int }_{-L}^L{\varphi }_h{\psi }_h\mathrm{d}x\forall {\vartheta }_h,{\psi }_h\in V_h.$$

(16)

We will show that problem (15) reduces to a system of linear equations. To achieve this, we assume that $M=\dim V_h<\infty$ and let $\{e_1,\dots ,e_M\}$ be a basis of $V_h$. Then, for each $n\in \{0,1,\dots ,N\}$, there exist ${\alpha }_1^{n+1},\dots ,{\alpha }_M^{n+1}\in \mathbb{R}$ such that

$${\vartheta }_h^{n+1}=\sum _{j=1}^M{\alpha }_j^{n+1}{e}_j.$$

(17)

We therefore write (15) as follows: For each $n\in \{0,1,\dots ,N\}$, find ${\alpha }_1^{n+1},\dots ,{\alpha }_M^{n+1}\in \mathbb{R}$ such that

$$\sum _{j=1}^M{\alpha }_j^{n+1}\{B(e_j,e_i)+\mathsf{\Delta }{t}^{n}A(e_j,e_i)\}=\sum _{j=1}^M{\alpha }_j^{n}B(e_j,e_i)\forall i=1,\dots ,M.$$

(18)

In this way, if we set $a_{ij}:=A(e_j,e_i)$ and $b_{ij}:=B(e_j,e_i)$, so that

$${\alpha }^{n+1}:=\left({\alpha }_j^{n+1}\right)\in {\mathbb{R}}^M,\mathrm{A}:=\left(a_{ij}\right)\in {\mathbb{R}}^{M\times M},\mathrm{B}:=\left(b_{ij}\right)\in {\mathbb{R}}^{M\times M},$$

the matrix form of (18) reads the following: For each $n\in \{0,1,\dots ,N\}$, find ${\alpha }^{n+1}\in {\mathbb{R}}^M$, which satisfies

$$(\mathrm{B}+\mathsf{\Delta }{t}^n\mathrm{A}){\alpha }^{n+1}=\mathrm{B}{\alpha }^n,$$

(19)

where the iteration is initialized with ${\alpha }^0\in {\mathbb{R}}^M$ obtained from (14). The following result establishes the unique solvability of the linear system (19).

Theorem 1.

The matrix $\mathrm{B}+\mathsf{\Delta }{t}^n\mathrm{A}$ is symmetric and positive definite, and therefore invertible.

Proof. 

The symmetry property follows directly from the definition of the bilinear forms A and B. Next, given $\beta :=\left({\beta }_j\right)\in {\mathbb{R}}^M$, we set

$${\psi }_h=\sum _{j=1}^M{\beta }_j{e}_j.$$

(20)

Proceeding analogously to ([36], Chapter 4), that is, using the $H_0^1\left(I\right)$-ellipticity of the form A with ellipticity constant $C>0$, depending on the diffusion coefficient D and the constant from Poincaré’s inequality (cf. (5)), we obtain

$$\begin{array}{cc}\hfill {\beta }^{\mathrm{T}}(\mathrm{B}+\mathsf{\Delta }{t}^n\mathrm{A})\beta & =\sum _{i,j=1}^M(b_{ij}+\mathsf{\Delta }{t}^n{a}_{ij}){\beta }_i{\beta }_j=B({\psi }_h,{\psi }_h)+\mathsf{\Delta }{t}^{n}A({\psi }_h,{\psi }_h)\hfill \\ \hfill & \ge \parallel {\psi }_h{\parallel }_{0,I}^2+\mathsf{\Delta }{t}^{n}C\parallel {\psi }_h{\parallel }_{1,I}^2\ge \mathsf{\Delta }{t}^{n}C{\parallel {\psi }_h\parallel }_{1,I}^2.\hfill \end{array}$$

Since ${\beta }^{\mathrm{T}}(\mathrm{B}+\mathsf{\Delta }{t}^n\mathrm{A})\beta >0$ for all $\beta$ values that are different from the null vector, the result follows.    □

2.3.2. Specific FE Subspace

This section specifies the matrix structure of the linear system (19) for a particular choice of $V_h$. Let ${\left\{x_m\right\}}_{0\le m\le M+1}$ be a uniform partition of the interval $\overline{I}=[-L,L]$, with the meshsize being denoted by $h>0$. We set

$$V_h=\{{\psi }_h\in C\left(\overline{I}\right):{{\psi }_h\mid }_{[x_{j-1},x_j]}\in {\mathcal{P}}_1\left([x_{j-1},x_j]\right)\forall j=1,\cdots ,M+2\}\cap H_0^1\left(I\right),$$

(21)

where ${\mathcal{P}}_1\left(S\right)$ denotes the space of polynomials of degree $\le 1$ defined over an interval S.

The following definition can be found in ([37], Section 1.1.2).

Definition 1.

For each $i\in \{1,\dots ,M\}$, the hat functions $e_i\in V_h$ are defined as

$$e_i\left(x\right)=\left\{\begin{array}{cc}{\displaystyle \frac{x-x_{i-1}}{h}}& x\in [x_{i-1},x_i],\hfill \\ {\displaystyle \frac{x_{i+1}-x}{h}}& x\in [x_i,x_{i+1}],\hfill \\ 0& x\notin [x_{i-1},x_{i+1}].\hfill \end{array}\right.$$

An example of a hat function is shown in Figure 2. It follows that $e_i\left(x_j\right)={\delta }_{ij}$, where ${\delta }_{ij}$ is the Kronecker delta function. Moreover, $\{e_1,\dots ,e_M\}$ is a basis for $V_h$; hence, a function ${\psi }_h\in V_h$ is uniquely determined by the values ${\left\{{\psi }_h\left(x_i\right)\right\}}_{1\le i\le M}$, namely

$${\psi }_h=\sum _{i=1}^M{\psi }_h\left(x_i\right)e_i.$$

Note that $a_{ij}=b_{ij}=0$ when $|i-j|\ge 2$. Furthermore, we obtain, after some algebraic manipulations,

$$a_{ij}=\left\{\begin{array}{cc}{\displaystyle \frac{2D}{h}}& \mathrm{if}j=i,\hfill \\ {\displaystyle -\frac{D}{h}}& \mathrm{if}|j-i|=1.\hfill \end{array}\right.$$

Similarly,

$$b_{ij}=\left\{\begin{array}{cc}{\displaystyle \frac{2h}{3}}& \mathrm{if}j=i,\hfill \\ {\displaystyle \frac{h}{6}}& \mathrm{if}|j-i|=1.\hfill \end{array}\right.$$

Consequently, the global matrix of (19) is tridiagonal, which is highly desirable in situations where the dimension of $V_h$ is very large, as this allows for a reduction in the number of flop required to solve the system. In particular, a tridiagonal system can be solved using the Thomas algorithm, which requires significantly fewer flop than the method that directly computes the inverse of the matrix. For further details, we refer the reader to [38].

2.3.3. Numerical Treatment of Non-Homogeneous Boundary Conditions

Recall from (9) that the unknowns $\varphi$ and $\vartheta$ are related by the decomposition $\varphi =\vartheta +\lambda$, where $\lambda$ is defined such that $\lambda ={\varphi }_{\infty }$ at $x=\pm L$. This approach allows $\varphi$ to be recovered from $\vartheta$, even in the presence of non-homogeneous boundary conditions. Building on this idea, we approximate $\varphi$ using the decomposition

$${\varphi }_h^{n+1}={\vartheta }_h^{n+1}+{\lambda }_h,n\in \{0,1,\dots ,N\},$$

(22)

where ${\vartheta }_h^{n+1}\in V_h$ is given by (17), with $V_h$ being defined in (21), and ${\lambda }_h$ is a polynomial satisfying the boundary conditions ${\lambda }_h(-L)={\lambda }_h\left(L\right)={\varphi }_{\infty }$. Specifically, we define

$${\lambda }_h:={\varphi }_{\infty }(e_0+e_{M+1}),$$

(23)

where $e_0$ and $e_{M+1}$ are the hat functions associated with the endpoints of I.

Next, using the same notation as in Section 2.3.2, the fully discrete formulation of (11) reads as follows: for each $n\in \{0,1,\dots ,N\}$, find ${\vartheta }_h^{n+1}\in V_h$ such that ${\vartheta }_h^0={\varphi }_h^0$ (cf. (14)) and

$$B({\vartheta }_h^{n+1},{\psi }_h)+\mathsf{\Delta }{t}^{n}A({\vartheta }_h^{n+1},{\psi }_h)=B({\varphi }_h^n,{\psi }_h)-\mathsf{\Delta }{t}^{n}D{\int }_{-L}^L{\lambda }_h^{\prime }{\psi }_h^{\prime }\mathrm{d}x\forall {\psi }_h\in V_h.$$

(24)

It is clear that (24) coincides with (15) when ${\varphi }_{\infty }=0$.

Now let $\mu \in {\mathbb{R}}^M$ be the vector whose entries are all zero, except for the first and last ones, which are both set to ${\varphi }_{\infty }$. Then, for each $n\in \{0,1,\dots ,N\}$, the discretization (24) yields

$$(\mathrm{B}+\mathsf{\Delta }{t}^n\mathrm{A}){\alpha }^{n+1}=\mathrm{B}{\alpha }^n-\mathsf{\Delta }{t}^n\mathrm{A}\mu .$$

(25)

Here, the matrices $\mathrm{A}$ and $\mathrm{B}$ are defined as in Section 2.3.2, thus ensuring the unique solvability of the linear system (25), which follows directly from Theorem 1.

We end this section by noting that the discrete concentration with non-homogeneous boundary conditions can be computed from (22) using the solution of (25), which provides

$${\varphi }_h^{n+1}=\sum _{j=0}^{M+1}{\alpha }_j^{n+1}{e}_j,n\in \{0,\cdots ,N\},$$

with ${\varphi }_h^{n+1}(-L)={\varphi }_h^{n+1}\left(L\right)={\varphi }_{\infty }$, as required.

2.4. Variable Diffusion Coefficient

This section presents an FE method to solve the non-linear equation in (1), assuming ${\varphi }_{\infty }$ is time-dependent for greater generality. Although this model is more challenging to analyze with FE methods due to the concentration dependence of D, it still relies on the previous results.

The weak form of (1) resembles (11), with the form A now replaced by

$$A(\varphi ;\vartheta ,\psi ):={\int }_{-L}^{L}D\left(\varphi \left(t\right)\right){\vartheta }^{\prime }\left(t\right){\psi }^{\prime }\mathrm{d}x\forall \psi \in H_0^1\left(I\right).$$

Note here that $\vartheta \left(t\right):=\varphi \left(t\right)-\lambda \left(t\right)$ and $\lambda \left(t\right)\in H_0^1{\left(I\right)}^{\perp }$ with $\lambda \left(t\right)={\varphi }_{\infty }\left(t\right)$ at $x=\pm L$.

Proceeding analogously to Section 2.3, we discretize the weak formulation of (1) using an FE method in space and an implicit Euler method in time. To handle the non-linearities, we employ a Picard-type iteration. Specifically, we consider the following fully discrete scheme: for each $n\in \{0,1,\dots ,N\}$ and given ${\varphi }_h^n$, find ${\vartheta }_h^{n+1}\in V_h$ such that ${\vartheta }_h^0={\varphi }_h^0$ and

$$\begin{array}{cc}\hfill & B({\vartheta }_h^{n+1},{\psi }_h)+\mathsf{\Delta }{t}^{n}A({\varphi }_h^n;{\vartheta }_h^{n+1},{\psi }_h)\hfill \\ \hfill & =B({\vartheta }_h^n,{\psi }_h)-{\int }_{-L}^L{\lambda }_h^{n+1}{\psi }_h\mathrm{d}x-\mathsf{\Delta }{t}^n{\int }_{-L}^{L}D\left({\varphi }_h^n\right){\left({\lambda }_h^{n+1}\right)}^{\prime }{\psi }_h^{\prime }\mathrm{d}x\forall {\psi }_h\in V_h,\hfill \end{array}$$

(26)

where ${\varphi }_h^0$ is given by (14) and B by (16). Furthermore, the last two integrals in (26) arise from the treatment of the time-dependent boundary condition. In fact, proceeding as in Section 2.3.3, we obtain

$${\varphi }_h^{n+1}={\vartheta }_h^{n+1}+{\lambda }_h^{n+1},\mathrm{with}{\lambda }_h^{n+1}:={\varphi }_{\infty }^{n+1}(e_0+e_{M+1}).$$

(27)

Now, inspired by [32], we define

$$D\left(\varphi \right):=D_0{(1+\delta \varphi )}^k,$$

(28)

where $k\in \mathbb{N}$, $\delta \in \mathbb{R}$, and $D_0>0$ are experimentally determined constants. For the sake of simplicity, we assume that k is either 1 or 2. However, the analysis in this section can be easily extended to cases where $k>2$. We furthermore assume that D is bounded: There exists a constant $D_{\ast }>0$ such that, for all $\psi \in [0,1]$,

$$D\left(\psi \right)\ge D_{\ast }.$$

(29)

Remark 1.

As discussed in [39], a general form for the diffusion coefficient is $D\left(\varphi \right)=D_{0}F\left(\varphi \right)$, where $D_0=D\left(0\right)$ is a positive constant. The choice in (28), with $F\left(\varphi \right)={(1+\delta \varphi )}^k$, represents a specific instance of this general form. While other forms of F, such as exponential functions, have been explored in the literature (see, e.g., [40]), we restrict our focus to (28) for brevity.

Next, we specify the matrix form of the fully discrete scheme (26), subjected to the diffusion coefficient (28). Let $\widehat{\mathrm{A}}:=\left({\widehat{a}}_{ij}\right)\in {\mathbb{R}}^{M\times M}$ denote the matrix given by

$${\widehat{a}}_{ij}:=A({\varphi }_h^n;e_j,e_i)={\int }_{-L}^{L}D\left({\varphi }_h^n\right)e_i^{\prime }{e}_j^{\prime }\mathrm{d}x=D_0{\int }_{-L}^L({1+\delta \sum _{l=0}^{M+1}{\alpha }_l^n{e}_l)}^k{e}_i^{\prime }{e}_j^{\prime }\mathrm{d}x,$$

(30)

where $e_0,e_1,\dots ,e_{M+1}$ are the hat functions introduced in Section 2.3.2.

After some algebraic manipulations, we obtain

$$\widehat{\mathrm{A}}={\displaystyle \frac{D_0}{h}}\left[\begin{array}{ccccccc}{r}_1& q_1& 0& 0& \cdots & \cdots & 0\\ q_1& r_2& q_2& 0& \cdots & \cdots & 0\\ 0& \ddots & \ddots & \ddots & \ddots & & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& q_{M-3}& r_{M-2}& q_{M-2}& 0\\ 0& \cdots & \cdots & 0& q_{M-2}& r_{M-1}& q_{M-1}\\ 0& \cdots & \cdots & 0& 0& q_{M-1}& r_M\end{array}\right],$$

with matrix components depending on the value of k. Specifically, for $k=1$,

$$\begin{array}{cc}\hfill q_i& :=-1-{\displaystyle \frac{\delta }{2h}}({\alpha }_i^n+{\alpha }_{i+1}^n),\hfill \\ \hfill r_i& :=2+{\displaystyle \frac{\delta }{2h}}({\alpha }_{i-1}^n+2{\alpha }_i^n+{\alpha }_{i+1}^n),\hfill \end{array}$$

while for $k=2$,

$$\begin{array}{cc}\hfill q_i& :=-{\displaystyle \frac{1}{3}}{(1+\delta {\alpha }_i)}^2-{\displaystyle \frac{1}{3}}(1+\delta {\alpha }_i)(1+\delta {\alpha }_{i+1})-{\displaystyle \frac{1}{3}}{(1+\delta {\alpha }_{i+1})}^2,\hfill \\ \hfill r_i& :={\displaystyle \frac{1}{3}}{(1+\delta {\alpha }_{i-1})}^2+{\displaystyle \frac{2}{3}}{(1+\delta {\alpha }_i)}^2+{\displaystyle \frac{1}{3}}{(1+\delta {\alpha }_{i+1})}^2+{\displaystyle \frac{1}{3}}(1+\delta {\alpha }_i)(2+\delta {\alpha }_{i-1}+\delta {\alpha }_{i+1}).\hfill \end{array}$$

Finally, for each $n\in \{0,1,\cdots ,N\}$, we set ${\mu }^{n+1}\in {\mathbb{R}}^M$ as the vector whose components are all zero, except for the first and last ones, which are set to ${\varphi }_{\infty }^{n+1}$. Then, the matrix form of (26) reads as follows: for each $n\in \{0,1,\dots ,N\}$, find ${\alpha }^{n+1}:=\left({\alpha }_j^{n+1}\right)\in {\mathbb{R}}^M$ such that

$$(\mathrm{B}+\mathsf{\Delta }{t}^n\widehat{\mathrm{A}}){\alpha }^{n+1}=\mathrm{B}{\alpha }^n-(\mathrm{B}+\mathsf{\Delta }{t}^n\widehat{\mathrm{A}}){\mu }^{n+1},$$

(31)

where the matrix $\mathrm{B}$ is defined as in Section 2.3.2.

Remark 2.

The unique solvability of the fully discrete scheme (31) follows from the fact that the diffusion coefficient (28) satisfies the condition (29). In fact, the matrix $\mathrm{B}+\mathsf{\Delta }{t}^n\widehat{\mathrm{A}}$ is symmetric, and for any vector $\beta \in {\mathbb{R}}^M$ given by (20), we have

$${\beta }^T(\mathrm{B}+\mathsf{\Delta }{t}^n\widehat{\mathrm{A}})\beta \ge \mathsf{\Delta }{t}^{n}C{\parallel {\psi }_h\parallel }_{1,I}^2,$$

with constant $C>0$ depending on $D_{\ast }$ and the constant from Poincaré’s inequality. It then follows that the matrix $\mathrm{B}+\mathsf{\Delta }{t}^n\widehat{\mathrm{A}}$ is positive definite, and therefore invertible.

3. Results

3.1. Drug Release Experiments

Samples of NC/nPSi composites (0.5 × 0.5 ${\mathrm{cm}}^2$) with varying concentrations of microparticulate nPSi were synthesized following the protocol outlined by K. Garrido-Miranda et al. [12]. Each sample had different thicknesses based on the percentage of nPSi (m/m): NC (control, 0.0%) = 6.5 ± 1 $\mathsf{\mu }\mathrm{m}$, NC/nPSi-0.1% = 10.5 ± 2 $\mathsf{\mu }\mathrm{m}$, NC/nPSi-0.5% = 12.7 ± 3 $\mathsf{\mu }\mathrm{m}$, and NC/nPSi-1.0% = 29.5 ± 4 $\mathsf{\mu }\mathrm{m}$ (Figure 3). Subsequently, the samples were loaded with a concentrated MB solution (0.001 M, pH 7.0) for 15 min at 100 rpm. They were then rinsed with distilled water and dried at room temperature. Release profiles were conducted in vials filled with 3 ml of saline phosphate-buffered solution (PBS, pH 7 and 37 °C) at 100 rpm on a horizontal shaker (INB-2005 LN, Biotek, Winooski, VT, USA). The concentration of MB in the fluid was measured at specific time intervals using UV-Vis spectrophotometry (UV-1800 Shimadzu, Kyoto, Japan) at a wavelength of 671 nm [7]. The release profiles were obtained as the mean of triplicate experiments (see the Supplementary Materials section).

3.2. Model Prediction

In this section, we compare the drug release experiments with the numerical results obtained using the FE method described in Section 2.4. All simulations were implemented using Octave 8.4.0 [33].

In what follows, the diffusion coefficient D is defined as in (28), with the constant $D_0$ determined using the Nelder–Mead optimization method [41], which minimizes the error

$$e_{ana}:=\sqrt{\sum _{n=0}^{N+1}{\left|u^n-w^n\left(D_0\right)\right|}^2},$$

where $u^n$ represents the experimental release profile at time $t^n$, and $w^n\left(D_0\right)$ denotes the corresponding value obtained from the analytical release profile in (3), computed as

$$w^n\left(D_0\right)=1-{\displaystyle \frac{8}{{\pi }^2}}\sum _{k=0}^p{\displaystyle \frac{1}{{(2k+1)}^2}}exp\left[-{\displaystyle \frac{D_0{(2k+1)}^2{\pi }^2{t}^n}{4L^2}}\right].$$

Note here that p is a user-defined integer, which we set to $p=100$. The computed values of $D_0$ for the four samples detailed in Section 3.1 are listed in

Table 1. Computed values of $D_0$.

Sample	${\mathit{D}}_0$ (${\mathsf{\mu }\mathbf{m}}^2$/h)
NC	28.648
NC/nPSi-0.1%	107.640
NC/nPSi-0.5%	17.413
NC/nPSi-1.0%	51.523

. Once the parameter $D_0$ is determined, we proceed to complete the definition of the diffusion coefficient D by selecting an appropriate value for the scalar $\delta$, as will be discussed in more detail below. With D being fully defined, the release profile $M_t/M_{\infty }$ is approximated by solving Fick’s second law using the FE method. Specifically, the solution to the system (31), as obtained from Algorithm 1, provides the following approximation:

$${{\displaystyle \frac{M_t}{M_{\infty }}\mid }}_{t=t^{n+1}}\approx {\displaystyle \frac{{\int }_{-L}^L({\varphi }_h^{n+1}-{\varphi }_0)\mathrm{d}x}{{\int }_{-L}^L({\varphi }_{\infty }^{n+1}-{\varphi }_0)\mathrm{d}x}}.$$

(32)

Unless otherwise specified, we choose ${\varphi }_0=1$ and ${\varphi }_{\infty }^{n+1}=0$ for all $n\in \{0,1,\dots ,N+1\}$.

Algorithm 1: FE solution.

Input: N, M, $\delta$, $D_{\mathrm{ref}}$, ${\varphi }_0$, ${\varphi }_{\infty }$, L      Output: ${\varphi }_h^0,{\varphi }_h^1,\dots ,{\varphi }_h^{N+1}$ 1: $h\leftarrow 2L/(M+1)$ 2: ${\varphi }_h^0\leftarrow$$L^2$-projection of ${\varphi }_0$ into $V_h$ (cf. (14)) 3: for $n=0,1,\dots ,N$ do 4:     $\mathsf{\Delta }{t}^n\leftarrow t^{n+1}-t^n$ 5:     ${\alpha }_0^{n+1}\leftarrow {\varphi }_{\infty }^{n+1}$ 6:     ${\alpha }_{M+1}^{n+1}\leftarrow {\varphi }_{\infty }^{n+1}$ 7:     ${\alpha }_1^{n+1},\dots ,{\alpha }_M^{n+1}\leftarrow$ solution to the problem (31) 8:     ${\varphi }_h^{n+1}\leftarrow {\sum }_{j=0}^{M+1}{\alpha }_j^{n+1}{e}_j$ 9: end for

Selecting the optimal values for $\delta$ and k is essential for accurately defining the diffusion coefficient D. To determine these values, we focused on the samples in Table 1. For each sample, k was set to either 1 or 2, the meshsize h was computed as $2L/(M+1)$ with $M=31$, and $\delta$ was chosen from the interval $[-1,1]$. The goal was to minimize the numerical error

$$e_{num}:=\sqrt{{\sum }_{n=0}^{N+1}{\left|u^n-v^n\left(\delta \right)\right|}^2},$$

(33)

where the experimental release profile at time $t^n$ is denoted by $u^n$, while $v^n\left(\delta \right)$ represents the corresponding approximation from (32). Values of $\delta$ outside the interval $[-1,1]$ led to larger errors and were therefore discarded. In this context, the computed errors for $k=2$ are shown in

Table 2. Release profile of errors associated with the variable diffusion coefficient D for $k=2$.

NC
$\delta$	−1	−0.3	−0.2	0	0.2	0.3	1
$e_{num}$	1.636	0.320	0.262	0.166	0.104	0.092	0.200
$e_{ana}$	–	–	–	0.140	–	–	–
NC/nPSi-0.1%
$\delta$	−1	−0.3	−0.2	0	0.2	0.3	1
$e_{num}$	1.573	0.333	0.279	0.194	0.145	0.135	0.193
$e_{ana}$	–	–	–	0.178	–	–	–
NC/nPSi-0.5%
$\delta$	−1	−0.3	−0.2	0	0.2	0.3	1
$e_{num}$	1.344	0.314	0.264	0.185	0.141	0.134	0.263
$e_{ana}$	–	–	–	0.171	–	–	–
NC/nPSi-1.0%
$\delta$	−1	−0.3	−0.2	0	0.2	0.3	1
$e_{num}$	1.219	0.306	0.257	0.191	0.176	0.185	0.351
$e_{ana}$	–	–	–	0.146	–	–	–

, indicating that the optimal value of $\delta$ is $0.3$ for the first three samples and $0.2$ for the last one. Thus, condition (29) is satisfied with $D_{\ast }=D_0$, ensuring the unique solvability of the fully discrete scheme (26), as stated in Remark 2. Similar results for $k=1$ were obtained and are omitted for brevity.

Figure 4 shows the experimental release data, their numerical approximations using the FE method, and the analytical solution from (3). The FE solution provided the best fit for the experimental data in the first three samples. For the last sample, where the concentration of nPSi was higher, the analytical release profile was 20.55% more accurate, as indicated by the errors in Table 2. It follows that D strongly depends on the concentration when the percentage of nPSi is below 1.0%. Additionally, the simulations in Figure 5 reveal that the diffusion rate decreases at higher nPSi percentages, as reflected by the slower reduction in concentration across the thickness, particularly in the last sample. These results align with the experimental findings of K. Garrido-Miranda et al. [12], where incorporating nPSi into the material enhances control over the release of MB.

4. Discussion

The FE method presented in Section 2.4 provides a general framework for accurately describing the drug delivery system under consideration. In fact, the variable diffusion coefficient defined in (28) showed strong agreement with data from K. Garrido-Miranda et al. [12], particularly for samples with nPSi percentages below 1.0%. A weaker dependency was observed at an nPSi concentration of exactly 1.0%. This may be linked to the increased thickness associated with higher nPSi levels, which was presumed negligible in both the numerical method and the analytical solution.

In our simulations, we assumed constant initial and boundary conditions. However, a more realistic assumption would involve time-dependent boundary data, which may occur when the solution is not well stirred [42]. Notably, the fully discrete scheme (26) remains effective in simulating this scenario, as illustrated in Figure 6 for a single sample, where the boundary data were derived from the empirical relation ${\varphi }_{\infty }\left(t\right)=1-M_t/M_{\infty }$. These results were contrasted with the analytical solution from (3), showing that the numerical solution with time-dependent ${\varphi }_{\infty }$ performs better when $M_t/M_{\infty }<0.85$. This conclusion is further supported by the error computed from (33) using $\delta =0.45$ and $k=2$, yielding a value of $0.149$, which represents a $14.77\%$ improvement compared to the corresponding error reported in Table 2 for the analytical solution with constant boundary conditions.

Building on these results, a more comprehensive approach to designing drug delivery systems that undergo significant geometric changes and time-dependent boundary conditions should incorporate space–time FE methods (see, e.g., [43,44]). Future work will explore this approach by incorporating additional parameters into the general form of diffusion coefficient presented in Remark 1, with the aim of capturing effects such as evaporation within the composite.

5. Conclusions

This study employed an FE method to model diffusion and controlled drug release from NC/nPSi composites based on Fick’s second law with variable diffusivity. The FE simulations, supported by experimental validation, demonstrated that increasing the nPSi concentration in the NC matrix enhances control over the release MB. This finding is particularly relevant for pharmaceutical applications, where controlled drug delivery is critical to minimize adverse effects on tissues. Furthermore, this study highlights the importance of considering geometric changes in the composite matrix, especially at high nPSi concentrations, which significantly influence drug release behavior. Addressing these complexities is a priority for future work to extend the applicability of the proposed model to a wider range of material configurations and drug delivery scenarios. The results not only validate the effectiveness of the FE method for modeling diffusion in composite materials but also provide insights into optimizing the design of NC/nPSi composites for controlled drug release. This work sets a foundation for further research into more complex geometries and dynamic environmental conditions, reinforcing its potential contributions to the development of advanced drug delivery systems.