Modeling the Kinetics, Curing Depth, and Efficacy of Radical-Mediated Photopolymerization: The Role of Oxygen Inhibition, Viscosity, and Dynamic Light Intensity

Introduction

Photoinitiated (photosensitized) polymerization and cross-linking provide advantageous means over thermal-initiated polymerization, including fast and controllable reaction rates, and spatial and temporal control over the formation of the material, without the need for high temperatures or harsh conditions (Fouassier, 1995; Odian, 2006). Tissue engineering using scaffold-based procedures for chemical modification of polymers has been reported to improve its mechanical properties by cross-linking or polymerization with ultraviolet (UV) or visible light to produce gels or high-molecular-weight polymers (Chen and Shi, 2014; Kotisch et al., 2017). The progress of light-responsive smart nanomaterials was recently review by Yang et al. (2018).

Industrial applications include developing materials for applications such as thin films, coatings, printing, graphic work, dentistry, contact lenses, and electronics. It is a noncontact, low-energy, and rapid process with capabilities of spatially specifying the reaction via photomasks (photolithography) (Cabral et al., 2004; O'Brien and Bowman, 2006; Cramer et al., 2008; Dendukuri et al., 2008; Alvankarian and Majlis, 2015; Wohlers and Caffrey, 2016; Chen et al., 2017; Wu et al., 2018). The polymerization rate is inhibited by air due to oxygen inhibition, which scavenges the radical species needed for cross-linking initialization. Diffusion of oxygen from a high-concentration zone into a prepolymer resin during UV curing requires an additional amount of photoinitiator and UV energy to consume the dissolved and diffused oxygen (O'Brien and Bowman, 2006; Cramer et al., 2008; Dendukuri et al., 2008; Alvankarian and Majlis, 2015).

Cabral et al. (2004) investigated the concept of frontal photopolymerization, in which the position of the frontal border can be controlled by adjusting UV power and the available atmospheric air. The utilization of microfabrication to reduce the deposition steps and to obtain a monolithic product was reported by Alvankarian and Majlis (2015). Chen et al. (2017) proposed a kinetic model for pillar growth that includes free-radical generation and oxygen inhibition in thick films of photoinitiated media and have demonstrated control over the pillar spacing and pillar height with the irradiation intensity, film thickness, and the size and spacing of the optical beams. In microfabrication system, the formation of radical decreases over depth due to the reduction in light intensity and photosensitizer (PS) concentration and increase in oxygen inhibition. Oxygen diffuses into the film and consumes radicals. The balance of radical production and oxygen inhibition gives rise to the inhibition zone, where the polymerization is completely suppressed (Wu et al., 2018).

The kinetics of polymerization has been extensively studied (Cabral et al., 2004; O'Brien and Bowman, 2006; Cramer et al., 2008; Dendukuri et al., 2008; Alvankarian and Majlis, 2015; Semchishen et al., 2015; Lin and Wang, 2016; Chen et al., 2017; Lin and Cheng, 2017; Lin, 2018, 2019; Wu et al., 2018; Lin et al., 2019). However, most of the previous models (Cabral et al., 2004; O'Brien and Bowman, 2006; Cramer et al., 2008; Dendukuri et al., 2008; Alvankarian and Majlis, 2015; Chen et al., 2017; Wu et al., 2018) are based on oversimplified assumptions of constant PS concentration (without depletion), and the light intensity in the polymer following a conventional Beer–Lambert law (BLL), with neglected PS depletion, which is only valid for optically thin polymers (Lin and Cheng, 2017).

This study will present a more general kinetic model with coupled equations for the PS and oxygen concentration and monomer, which will be numerically solved for various conditions. Analytic formulas for the efficacy and curing depth will be derived and utilized to analyze the measured pillar height in microfabrication reported by Chen et al. (2017). Various scaling laws for the efficacy and curing depth, in both transient and steady state, will be derived, for the first time. The importance of PS depletion and dynamic light intensity in an optically thick polymer, based on a revised BLL, will be shown numerically (Lin and Cheng, 2017). Finally, various strategy for improved efficacy and curing depth, by reducing oxygen inhibition effect, will be discussed.

Methods

Kinetic Equations

Photopolymerization in general includes radical-mediated, cationic, and anionic catalyzed, and atom transfer radical polymerization. In a radical-mediated photopolymerization, the monomer is converted to polymer after the light irradiation of the photosensitizer (PS) or photoinitiator. The UV (or visible) light-produced triplet state (T^*) is coupled to the substrate monomer [A] and the oxygen to produce radicals. Each radical becomes the center of origin of a polymer chain. The chain growth of a polymer radical with m-links stops as a result of chain termination reactions. Kinetic equations of an m-component radical photopolymerization process (with a triplet excited state as the catalyst) may be described as follows: (Semchishen et al., 2015; Lin, 2019)

$\begin{array}{ll}\frac{\partial T(z,t)}{\partial t}=bIC-(k_5+k_3[{\text{O}}_2]+k_7[\text{A}])\text{T}& (1)\end{array}$

$\begin{array}{ll}\frac{\partial [{\text{O}}_2]}{\partial t}=P-k_{3}T[{\text{O}}_2]+-{\displaystyle {\sum }_{m=1}^{\infty }{k}_m{R}_m[{\text{O}}_2]}& (2)\end{array}$

$\begin{array}{llll}\frac{\partial R_1}{\partial t}& =& k_8\left[\text{A}\right]\text{T}-{\text{k}}^{\prime }{\text{R}}_1\left[\text{A}\right]-{\text{R}}_1{\sum }_{\text{m}=1}^{\infty }{\text{k}}_{\text{m}}{\text{R}}_{\text{m}}& (3)\end{array}$

$\begin{array}{llll}\frac{\partial R_{n+1}}{\partial t}& =& k^{\prime }(R_n-R_{n+1})\left[\text{A}\right]-R_{n+1}{\sum }_{m=1}^{\infty }{k}_m{R}_m& (4)\end{array}$

$\begin{array}{llll}\frac{\partial \left[\text{A}\right]}{\partial t}& =& -k_7\left[\text{A}\right]T-{\sum }_{m=1}^{\infty }{k}_m{R}_m\left[\text{A}\right]& (5)\end{array}$

As shown in Figure 1, for a one-component monomer A, with three-radical system consisting of two PS radicals (R′ and R) and one singlet oxygen, the pathways are described as follows. The ground-state PS molecules are excited by the UV light to its singlet excited state (S₁), which could be relaxed to its ground state or to a triplet excited state (T^*). In type-I process, T^* could interact directly with the substrate [A] and produces the first-radical (R′), which could produce (by chain reaction) a second radical (R) which could interact with the ground state oxygen or the first radical. For type-II process, (T₃) interacts with [O₂] to form oxygen singlet [¹O₂]. These reactive radicals, R, R′, and [²O₁], could be relaxed to its ground states or interacts with the substrate [A]. This article will focus on the one-monomer system, and the two-monomer and two-initiator systems will be shown elsewhere (Chen et al., 2019b; Lin, 2019).

FIGURE 1

Figure 1. Schematics of photochemical pathways. Pathway of one-monomer system, for radical-mediated pathway-1 and oxygen-mediated pathway-2; where P is the ground-state photosensitizer, having an excited and triplet state PS* and T*, which interacts with the substrate A to form radicals R′ and R.

For a one-monomer system, using the short-hand notations, the concentration of various components: C (z, t) and T for the PS ground and triplet state, respectively; [O₂] and X for the oxygen ground and singlet state; R and R′ for the first and intermediate free radical; and [A] for the available extracellular matrix substrate; and the kinetic equations for a three-radical system becomes

$\begin{array}{ll}\frac{\partial C(z,t)}{\partial t}=-bI\left(z,t\right)C-k_{11}XC+(k_5+k_3[{\text{O}}_2])T+R_E& (6)\end{array}$

$\begin{array}{ll}\frac{\partial T}{\partial t}=bI\left(z,t\right)C-(k_5+k_3[{\text{O}}_2]+k_7[\text{A}])T& (7)\end{array}$

$\begin{array}{llll}\frac{\partial R^{\prime }}{\partial t}& =& 2k_{7}T\left[\text{A}\right]-{\text{k}}_{12}\text{R}{\text{R}}^{\prime }-{2\text{k}}_{\text{T}}{{\text{R}}^{\prime }}^2& (8)\end{array}$

$\begin{array}{l}\frac{\partial R}{\partial t}=2k_T{R^{\prime }}^2-kR\left[\text{A}\right]-k_{12}RR-2k_T{R}^2\\ -k^{″}R[{\text{O}}_2]& (9)\end{array}$

$\begin{array}{ll}\frac{\partial X}{\partial t}=k_3[{\text{O}}_2]T-(k_6+k_{11}C+k_8[\text{A}])X& (10)\end{array}$

$\begin{array}{ll}\frac{\partial [{\text{O}}_2]}{\partial t}=P-k_3[{\text{O}}_2]T+k_{6}X-k{R}^{″}[{\text{O}}_2]& (11)\end{array}$

$\begin{array}{ll}\frac{\partial [\text{A}]}{\partial t}=-(k_{7}T+k_{8}X+k^{\prime }R)[\text{A}]& (12)\end{array}$

where b = 83.6a′qw, w being the UV light wavelength (in cm) and light intensity I(z,t) in mW/cm²; q is the quantum yield of the PS triplet state; Equation (11) also includes an oxygen source term given by Semchishen et al. (2015), P = (1 – [O₂]/O₀)P′, with a maximum rate constant P′, with a maximum rate constant P₀ and initial oxygen concentration Y₀. This term may be also given by the oxygen diffusion P = $D_0{\nabla }^2\left[O\right].$ In Equation (1), the regeneration term is given by R_E = k″R[O₂] + 2k_T(R² + R^′2),

All the reaction rate constants are defined by the associated coupling terms, as shown by Figure 2. For examples, in Equation (9), k₇ is for the coupling rate of [A] and T, which has a ground-state relaxation rate k₅; in Equation (11), k′ is for [A] and R, which is coupled with R′ by k₁₂ and a bimolecular termination rate of k_T; in Equation (12), k′ is for the reaction of [A] and R, which is coupled with oxygen by k″ and a bimolecular termination rate k_T. Figure 2 also shows the rate constants for type-II process: k₃ for singlet-oxygen (X) production; k₈ for cross-link of X with [A] and having a relaxation rate of k₆; and a reduction rate k₁₁, by coupling to C.

FIGURE 2

Figure 2. The ground-state photosensitizer (PS) (C) is excited by a UV light to its singlet excited state (C₁), which could be relaxed to its ground state or to a triplet excited state (T*). In type-I process, T* could interact directly with the monomer [A] to generate free radical (R′) by recombination. The radical R could interact with [A] for cross-linking, or oxygen [O₂], or terminated by coupling with R′, or bimolecular recombination (2R²). For type-II process, T* interacts with [O₂] to form oxygen singlet [¹O₂] (X), which could relax to its oxygen [O₂], or interacts with [A] for cross-linking, or coupling with C. All rate constants are shown in reds next to arrows.

The dynamic light intensity is given by (Lin and Wang, 2016; Lin and Cheng, 2017)

$\begin{array}{llll}\frac{\partial I(z,t)}{\partial z}& =& -A^{\prime }(z,t)I(z,t)& (13)\end{array}$

$\begin{array}{llll}{A}^{\prime }(z,t)& =& 2.3\left[(a^{\prime }-b^{\prime })C(z,t)+b^{\prime }{C}_0+Q\right]& (14)\end{array}$

where a′ and b′ are the extinction coefficients of PS and the photolysis product, respectively; Q is the absorption coefficient of the monomer at the UV wavelength. Most previous modeling (Cramer et al., 2008; O'Brien and Bowman, 2006; Dendukuri et al., 2008; Alvankarian and Majlis, 2015; Chen et al., 2017; Wu et al., 2018) assumed a constant C (z, t) in Equation (14). Our analytic formulas in this article will use a time average of A (z, t) to count for the dynamic of light intensity due to PS depletion.

The kinetic Equations (6)–(12) may be numerically calculated to find the conversion efficacy, which, however, is too complex for us to analyze the roles of each of the parameters. For comprehensive modeling, we will use the so-called quasi-steady-state assumption (Lin, 2018; Lin et al., 2019). The life time of the singlet and triplet states of photosensitizer, the radicals (R and R′), and the singlet oxygen (X) are very short (nanosecond to microsecond time scale) since they either decay or react with cellular matrix immediately after they are created. Thus, one may set dT/dt = dR/dt = dR′/dt = dX/dt = 0, which give the quasi-steady-state solutions: T = bIgC, X = bIg′ [O₂]/k₈, g = k₈/(k₃[O₂] + k₈[A] + k₅); g′ = k₃/(k₆ + k₁₁C + k₇[A]);_. k = (k₁₁/k₈). The kinetic Equations (8)–(12) become,

$\begin{array}{ll}\frac{\partial C(z,t)}{\partial t}=-bI(g[\text{A}]+kg{g}^{\prime }C[{\text{O}}_2])C+R_E& (15)\end{array}$

$\begin{array}{ll}\frac{\partial [{\text{O}}_2]}{\partial t}=-(bICg+k^{″}R)[{\text{O}}_2]+P& (16)\end{array}$

$\begin{array}{ll}\frac{\partial [\text{A}]}{\partial t}=-{\left[bIgC\left(1+g^{\prime }\left[{\text{O}}_2\right]\right)+k^{\prime }R\right]}_{}\left[\text{A}\right]& (17)\end{array}$

where the PS regeneration term R_E = k₉R[O₂] + 2k_TR². The radicals, R and R′, are given by the solution of the following steady-state of Equations (8) and (9):

$\begin{array}{ll}2bIgC\left[\text{A}\right]-k_{12}R{R}^{\prime }-{2k}_T{R^{\prime }}^2=0& (18)\end{array}$

$\begin{array}{ll}2k_T{R^{\prime }}^2-k^{\prime }R\left[\text{A}\right]-k_{12}R{R}^{\prime }-2k_T{R}^2-k^{″}R[{\text{O}}_2]=0& (19)\end{array}$

Numerical solutions are required for R and R′. However, analytic formulas are available for the case that coupling of R′ and R is weaker than the bimolecular recombination, i.e., 2k_TR^′2 >> k₁₂RR′ in Equation (18), we obtain R′ = (bIgC[A]/k_T)^0.5. Substituting this steady-state R′ into Equation (19), we obtain

$\begin{array}{ll}R\left(z,t\right)=\left(\frac{1}{4k_T}\right){[-G+\sqrt{G^2+16k_{T}B(z,t)}]}_{}& (20)\end{array}$

where B = bICg[A], and G = k″ [O₂] + 2k₁₂R′ + k′[A]. For G² << 8k_TB, Equation (20) is further reduced to R = (B/k_T)^0.5 – k″ [O₂](1 – B′)/(4k_T), with B′ = 0.5k″[O₂]/(8Bk_T)^0.5, which shows that R, and efficacy, are decreasing functions of k″ [O₂], referred as the oxygen inhibition effect.

Equations (15)–(17) may be numerically solved. However, it requires all the rate constants (k_j) to be given for specific materials. Without knowing these parameters, we will further simplify Equations (15)–(17) as follows. In most situations, the monomer concentration ([A]) is much larger than the oxygen, which is also depleted faster than PS (or C), or k₈[A] >> (k₃[O₂] + k₅), such that g = 1/[A], k″[O₂]C/A² << 1, and bIC/A << k′R. These conditions were also reported by Dendukuri et al. (2008), Alvankarian and Majlis (2015), and Chen et al. (2017) in their modeling. Therefore, the depletion of C of Equation (15) is mainly due to type-I mechanism, whereas oxygen depletion is mainly due to the coupling term of k″R[O₂] in Equation (16). Equations (15)–(17) reduce to the following

$\begin{array}{llll}\frac{\partial C}{\partial t}& =& -bI(z,t)C(z,t)& (21)\end{array}$

$\begin{array}{ll}\frac{\partial [\text{O}]}{\partial t}=-k^{″}R\left(z,t\right)[{\text{O}}_2]+P& (22)\end{array}$

$\begin{array}{llll}\frac{\partial \left[A\right]}{\partial t}& =& {-k^{\prime }R(z,t)\left[\text{A}\right]}_{}& (23)\end{array}$

The above simplified kinetic equations reduce to that of Dendukuri et al. (2008) for the special situations (or assumptions): the PS concentration, C(z,t), is a constant, or Equation (21) dC/dt = 0 (or bIt << 1, for small dose), and light intensity is also a constant, or Equation (13), dI/dz = 0, which, however, is valid only for a short exposure time, or an optically thin polymer (with Az << 0.1). Wu et al. (2018) included the light intensity reduction in their thick polymer system. However, they have ignored the PS depletion, or dC/dt = 0. Furthermore, the rate constants, k′ and k″, in general, are reduction function of the monomer conversion efficacy due to the viscosity effects, which were also ignored by Dendukuri et al. (2008) and Alvankarian and Majlis (2015). They also used a simplified format of Equation (16) for the oxygen inhibition effect. Our Equations (15)–(17), to be numerically solved later, are more accurate than previous model (Dendukuri et al., 2008; Alvankarian and Majlis, 2015; Semchishen et al., 2015; Chen et al., 2017; Wu et al., 2018), which used our simplified Equations (22) and (23), with I(z,t) and C(z,t) are constants.

Accurate solutions of Equations (21)–(23) require numerical simulations (to be shown later). For analytic formulas, we will use approximated analytic formulas for the light intensity and the PS concentration, such that we do not need to solve for Equation (13), and the expressive closed forms of I(z,t) and C(z,t) allow us to solve for the first- and second-order solutions of [O₂], [A], and the conversion efficacy.

Analytic Formulas for Efficacy

The monomer conversion efficacy for a bimolecular termination process is given by C_EFF = 1 – [A]/[A]₀ = 1 – exp(–S), with [A₀] being the initial monomer concentration, and the S-function is given by the time integral of k′R. Without the oxygen diffusion (with P′ = 0), the first-order solution of Equation (22), with R = (0.5B/k_T)^0.5 is given by, [O]⁽¹⁾ = Y₀ exp(–S), Y₀ is the initial oxygen, which gives the second-order solution of Equation (6), R(2) = B/$k_{\text{T}}^{0.5}$ – k″ [O₂](1 – B′)/(4k_T), with B = bIC, B′ = 0.5k″ [O₂]/(Bk_T)^0.5, for B << 1. R′ (2) may be used to find the second-order solution of oxygen concentration (for P′ = 0) given by

$\begin{array}{ll}[{\text{O}}_2]=Y_0\exp (-S^{\prime })& (24)\end{array}$

$\begin{array}{llll}{S}^{\prime }& =& k^{\prime }{\displaystyle {\int }_0^t}{R}^{\prime }(2)dt& (25)\end{array}$

To include the oxygen diffusion effect, for small time when oxygen is present, the approximate solution of Equation (22) is given by (Dendukuri et al., 2008)

$\begin{array}{ll}[{\text{O}}_2]=Y_0-4k^{\prime }Bt/(\text{π}k)& (26)\end{array}$

The induction time (T_ID) is defined by when [O₂] = 0. We obtain

$\begin{array}{ll}{T}_{ID}=\pi k{Y}_0/(4k^{\prime }B)& (27)\end{array}$

During the induction time (with t < T_ID), the solution of Equation (23) is given by

$\begin{array}{llll}\left[\text{A}\right]& =& A_{0}exp(-S^{\prime })& (28)\end{array}$

$\begin{array}{llll}{S}^{\prime }& =& k^{\prime }{\displaystyle {\int }_{T_{ID}}^t}(B/G^{\prime })d{t}^{\prime }& (29)\end{array}$

where G′ = D₀Y₀(1 – 4k′Bt/(πk).

For large time when oxygen is largely depleted, the radical concentration is given by R = (B/k_T)^0.5 – K₁₂[O₂], with K₁₂ = k′(1 – B′)/(4k_T), with [O₂] given by Equation (24), which gives the solution of Equation (23),

$\begin{array}{llll}\left[\text{A}\right]& =& A_{0}exp(-S)& (30)\end{array}$

$\begin{array}{ll}S=k^{\prime }{\displaystyle {\int }_0^t(\sqrt{0.5B(z,t)/k_T}}-K_{12}[{\text{O}}_2])d{t}^{\prime }& (31)\end{array}$

where [O₂] may be analytically given by Equations (24) or (26). Equation (30) was also presented by Alvankarian and Majlis (2015). However, they assumed a constant PS concentration, or B is time independent. When the diffusion and consumption of oxygen inside the photo resin, or d[O₂]/dt = 0, the closed form solution is available. For the case that oxygen diffusion is stopped, or d[O₂]/dz = 0, the closed form solution of oxygen may be used to find the closed form solution of Equation (22) and the conversion efficacy. We note that, during the induction time (with t < T_ID), S′ of Equation (29) is much smaller than S of Equation (31); therefore, the conversion efficacy is a fast-rising function of time after t > T_ID. We will show the numerical results later.

The Dynamic Light Intensity

Solving Equations (13) and (14) for the light intensity, I(z,t) and PS concentration C(z,t) concentration, we may numerically find S′ and then the conversion efficacy. We may further derive the analytic form of conversion efficacy which requires a closed form of I(z,t) and C(z,t) as follows. Using our previously developed approximated analytic formulas (Lin and Wang, 2016; Lin and Cheng, 2017)

$\begin{array}{llll}I(z,t)& =& I_{0}exp\left[-A^{\prime }z\right]& (32)\end{array}$

$\begin{array}{llll}{C}^{\prime }(z,t)& =& C_{0}exp\left[-B^{\prime }t\right]& (33)\end{array}$

$\begin{array}{llll}{A}^{\prime }(z,t)& =& 2.3(a^{\prime }{C}_0+Q)-A_{1}t& (34)\end{array}$

where B′ = bI₀exp(–A″z), A₁ = 2.3(a′ – b′)C₀I₀bz, with A″ being the averaged absorption given by A″ = 0.5 × 2.3(a′ + b′) + 2.3Q. We note that the –A₁t term represents the decrease in A′ or increase in light intensity due to PS depletion. Using Equations (32)–(34), we obtain analytic solution of Equation (30), when K₁₂[O₂] is ignored in Equation (31) for a type-I predominated mechanism, as follows.

$\begin{array}{llll}S& =& KG(z,t)\sqrt{0.5bX{I}_0{C}_0}& (35)\end{array}$

$\begin{array}{llll}G(z,t)& =& \left[1-exp[-B^{″}{t}^{″})\right]/B^{″}& (36)\end{array}$

where K = k′/k_T^0.5, t″ = t – T_ID, B″ = 0.5(B′ – 0.5A₁t), X = exp[–A₂z], with A₂ = 2.3(a′C₀ + Q). We note that Equation (34) defines the dynamic feature of the light intensity which is an increasing function of time due to the depletion of the PS concentration. It also provides the non-linear dynamic dependence of A′(z,t), given by A₁t, which is important for optical-thick polymer and for high light dose. Greater detail analysis will be given later. The above analytic formulas provide useful information to analyze and predict critical roles of each of the influencing factors without numerically solving the coupled equations. Equation (35) will also be used to analyze the numerical results and the experimental data later.

The Effects of Viscosity

To include the viscosity effects, the free-volume is reduced when cross-link efficacy increases. As proposed by Wu et al. (2018), the rate parameters k′ and k_T, are decreasing function of the efficacy (C_EFF) given by (in CGS units, cm³/mM/s): k′ = 1,865/(1 + 2E′ × 10⁻⁹), k_T = 10⁷/(2.3 + E′) + k′(1 – C_EFF); with E′ = exp(vC_EFF), where v is a constant ranging 25–40 defining the viscosity strength. High viscosity effect leads to a lower conversion efficacy. To include the viscosity effect, the effective rate constant, K = k′/k_T^0.5, of Equation (14) may be revised to a smaller value as: 1 – m [1 – exp(–S)], with S is the first-order solution (with no revision), where m = 0.1–0.3 is a fit parameter to measured data. The revised K leads to a smaller efficacy due to the less free volume of cross-link resulted by the increase in viscosity when efficacy increases. Numerical results will be shown later.

Curing Depth and Inhibition Zone

A curing depth (Z_C) is defined by when the conversion efficacy is higher than a critical value, C_EFF > C_T, or when S > S_T, with S_T = ln [1/(1 – C_T)]. Using Equation (35), and let S = S_T = 2 (or C_T = 0.86), Z_C is related to the cross-link time (T_C) by, for the case that K₁₂[O₂] is ignored in Equation (31),

$\begin{array}{ll}{T}_C=\left(\frac{1}{B^{″}}\right)\ln [2B^{″}/(K^{\prime }\sqrt{0.5b{X}^{\prime }{I}_0{C}_0})-1]& (37)\end{array}$

where K′ = k′/k_T^0.5, X′ = exp(–A″Z_C) defines the curing depth (Z_C). We plot the curve of T_C vs. Z_C, then rotate the axis to show curve of Z_C vs. T_C. The above curing depth (Z_C) is proportional to the pillar height measured by Chen et al. (2017) where their measured data of Figure 4 will be analyzed by Equation (37) later. For a given photo resin thickness of H, Chen et al. (2017) also defined an inhibition zone (Z_D) given by Z_D = H – Z_C. When the oxygen inhibition effect, or k[O₂] is included in Equation (31), the analytic formula of Equation (37) is not available and Z_C needs a numerical solution to be shown later.

Results and Discussion

In the following numerical solutions of Equations (21)–(23), we will also use the radical R given by Equation (20) and light intensity given by Equation (32). The roles of key parameters of b, k″, I₀, and the initial PS and oxygen concentrations, polymer thickness, and the viscosity effect will be analyzed.

Dynamic Concentration and Light Intensity

Figure 3 shows the dynamic profiles of the PS concentration and light intensity obtained by the numerical solution of Equations (21)–(23), using light intensity given by Equation (32). Depending on the coupling parameter b and I(z,t), as shown by Equation (32), the depletion of C(z,t) causes the increasing profile of I(z,t), defined by a ratio R_Z′ = I(z,t)/I(z,t = 0), which is also a decreasing function of the depth (z) per BLL. Figure 3 shows the ratio R_Z, at z = 150 μm, for b = 0.01–0.02, a′ = 458 (1/mM), b′ = Q = 0. Figure 2 shows that higher intensity results a faster PS depletion and hence larger light intensity increase. For thick polymers (>100 μm). Therefore, the assumption of time-independent light intensity and PS concentration is valid only for optically thin (with A′z < 0.2) polymers and under a small dose, i.e., under the condition that the time-dependent term A₁t of Equation (34) is neglected. For optically thick polymers under a larger dose, with bI₀tz > 0.2, the thin-polymer assumption (O'Brien and Bowman, 2006; Cramer et al., 2008; Dendukuri et al., 2008; Alvankarian and Majlis, 2015; Chen et al., 2017; Wu et al., 2018) will cause an error of 10–30%, depending on the value of A′ in Equation (34). Figure 3B shows the dynamic feature of I(z,t)/I₀ due to PS depletion for z > 0.

FIGURE 3

Figure 3. Dynamic profiles of the PS concentration (blue curves) and light intensity increase ratio (red curves) at various light intensity I₀ = (0.5, 2, 5, 10) mW/cm², at z = 0 (A) and 150 μm (B).

Efficacy Temporal Profiles

We will now show the numerical results of the efficacy by solving Equations (21)–(23) using light intensity given by Equation (32). We will also show the dynamic profiles of oxygen. Roles of key parameters of b, k″, and I₀ and the initial PS and oxygen concentrations, polymer thickness, and the viscosity effect will be analyzed.

As shown by Figure 4, higher light intensity (I₀) and coupling factor (b) cause faster oxygen depletion (and shorter induction time) leading to higher conversion efficacy, as also predicted by Equations (24)–(35). Figure 5 shows that higher PS initial concentration (C₀) leads to higher efficacy, as predicted by our analytic Equation (35).

FIGURE 4

Figure 4. Dynamic profiles. Efficacy and oxygen profiles for various light intensity I₀ = (0.5, 2, 5, 10) mW/cm², for curves (1, 2, 3, 4), for b = 0.02, shown by (A,C); and b = 0.04, shown by (B,D); for C₀ = 0.01 mM, and without viscosity effect (or v = 0), and without external oxygen supply (with P′ = 0).

FIGURE 5

Figure 5. Dynamic profiles. Same as Figure 4, but for (A) C₀ = 0.1% and (B) C₀ = 0.2%; for b = 0.02.

Figure 6 shows that larger initial oxygen concentration causes higher inhibition effects and leads to lower efficacy, as predicted by Equation (31). It also shows longer induction time for lower light intensity and/or higher oxygen initial concentration. The induction time is defined by when oxygen is depleted.

FIGURE 6

Figure 6. Dynamic profiles. Same as Figure 4, but for oxygen initial value of 0.001 mM shown by (A,C); and 0.002 mM shown by (B,D); for b = 0.02 and C₀ = 0.01 mM.

Figure 7 shows that larger viscosity factor (v) causes smaller free volume and thus smaller radical coupling rate (k′) and leads to a lower efficacy, as predicted by the effective rate constant, K = k′/k_T^0.5, in Equation (35). Figure 8 shows that, inside the polymer (with z = 150 μm), both light intensity and PS concentration are reduced; thus, conversion is less, although oxygen is depleted slower inside the polymer, as predicted by the z-dependent X = exp[–A₂z] in Equation (31).

FIGURE 7

Figure 7. Dynamic profiles. Same as Figure 4, but for viscosity factor v = 28 (A); and v = 40 (B); for b = 0.02 and C₀ = 0.01 mM.

FIGURE 8

Figure 8. Dynamic profiles. Same as Figure 3, but for z = 0 (on surface), shown by (A,C); and z = 150 μm, shown by (B,D); for b = 0.02 and C₀ = 0.01 mM.

Figure 9 shows that external oxygen supply cause increased radical inhibition and thus smaller conversion, as predicted by Equation (34). Figure 10 shows the efficacy profiles for various external oxygen supply (P′ > 0) that larger P′ causes higher oxygen inhibition effect and leads to smaller efficacy, as predicted by Equation (31).

FIGURE 9

Figure 9. Dynamic profiles. Same as Figure 4, but for without (P′ = 0) and with (P′ = 0.001) external oxygen supply, shown by (A,B), respectively; for b = 0.02 and C₀ = 0.01 mM.

FIGURE 10

Figure 10. Dynamic profiles. Same as Figure 9, but for P′ = (0, 0.0001, 0.001, 0.006), for curves (1, 2, 3, 4); for b = 0.02, I₀ = 5 mW/cm² and C₀ = 0.01 mM.

Efficacy Spatial Profiles

The efficacy spatial profiles will be shown based on the analytic formula of Equation (35), with ignored oxygen effects, in Equation (31). Figure 11 shows that efficacy is a decreasing function of z, but it is an increasing function of light exposure time, as predicted by Equation (35). Figure 12 shows that efficacy is an increasing function of light intensity. The efficacy is also an increasing function of PS concentration (C₀), as shown by Figure 13. When the oxygen inhibition effect is included, as shown by Equation (31), the above analytic formulas produced efficacy will be reduced. However, only numerical data are available. The modeling curves of Alvankarian and Majlis (2015) shown by their Figures 3F,G are not accurate because their Equation (5) for the conversion have over simplified the oxygen inhibition effect, shown by our Equation (17) by a reduction factor of (1 – [O₂]) C_EFF, and also assumed a constant PS concentration throughout the resin thickness. Therefore, their efficacy curves in their Figure 7G are overestimated, compared to our Figure 10, which includes the PS depletion.

FIGURE 11

Figure 11. Efficacy spatial profiles. Efficacy vs. depth (z) for t = (60, 90, 120, 180) s, for curves (1, 2, 3, 4); for a fixed light intensity I₀ = 15 mW/cm², C₀ = 0.1 mM, [O₀] = 0.001 mM, b = 0.002, A′ = 1,000 (1/cm), in the absence of oxygen, based on Equation (35).

FIGURE 12

Figure 12. Efficacy spatial profiles. Same as Figure 11, but for I₀ = (2, 8, 30, 100) mW/cm², for curves (1, 2, 3, 4); for t = 400 s.

FIGURE 13

Figure 13. Efficacy spatial profiles. Same as Figure 11, but for C₀ = (0.1, 0.2, 0.3, 0.4) mM, for curves (1, 2, 3, 4); for I₀ = 15 mW/cm², at t = 400 s.

FIGURE 14

Figure 14. Curing depth profile. Curing depth (Z_C) vs. light intensity (I₀) for t = (60, 300, 600) s, for curves (1, 2, 3); for C₀ = 0.01 mM, based on Equation (15) without oxygen inhibition.

Curing Depth Profile

In microfabrication system, formation of radical (R′') decreases over depth (z) due to the reduction in light intensity and PS concentration and increase in oxygen inhibition. When oxygen diffuses into the film and consumes radicals, the balance of radical production and oxygen inhibition gives rise to the inhibition zone, where the polymerization is completely suppressed, as presented by kinetic Equation (37).

Using the analytic formula, Equation (37), we are able to investigate the roles of PS and oxygen concentration, light intensity, and exposure time in curing depth. Figure 13 shows that curing depth (Z_C) is an increasing function of light intensity and exposure time. Figure 15 shows curing depth profiles vs. time and light dose, in which the intensity dependence is reversed in Figures 15A,B. This unexpected feature also shown in the efficacy profile, shown by Equation (35), that smaller light intensity achieves higher steady-state efficacy (for a fixed light dose). Figure 16 shows that curing depth is an increasing function of the PS concentration (C₀), as also predicted by Equation (35) that efficacy has a scaling law of $C_0^{0.5,}$ and so does Z_C. The effect of time-dependent absorption factor A′ shown in Equation (35) on the curing depth is shown in Figure 17 that assumption of a constant PS concentration, or A₁t = 0 in Equation (34), will underestimate the cross-link depth 1–12%, especially for high dose and optically thick polymer. Figure 17 shows the viscosity effect on the efficacy. The effective rate constant (K) of Equation (35) is revised to a smaller value as: 1 – m [1 – exp(–S)], with S is the first-order solution (with no revision), where m = 0.1–0.3 is a fit parameter to measured data. The revised K leads to a smaller efficacy due to the less free volume of cross-link resulted by the increase of viscosity when efficacy increases.

FIGURE 15

Figure 15. Curing depth profile. Curing depth for I₀ = (5, 10, 15, 20) mW/cm², for curves (1, 2, 3, 4); for (A) Z_C vs. time and (B) Z_C vs. light dose.

FIGURE 16

Figure 16. Curing depth profile. Curing depth vs. time, for C₀ = (0.1, 0.2, 0.3, 0.4) mM, for curves (1, 2, 3, 4); for I₀ = 10 mW/cm².

FIGURE 17

Figure 17. Curing depth profile. Curing depth vs. time, with (blue) and without (red) the correction of time-dependent absorption factor A′ shown in Equations (13) and (14) for C₀ = 0.1 mM, I₀ = 10 mW/cm².

It is important to note that the S-function is scaled by [C₀I₀t²]^0.5 for transient state (with B′t << 1), and scaled by [C₀/I₀]^0.5 for steady state. Therefore, the dose-dependence Bunsen Roscoe law, scaled by a light dose E₀ = I₀t, failed in a bimolecular termination system, with radical (R′) scaled by [C₀I₀ exp(–B′t)]^0.5, as shown by Equation (35). These features also apply to the curing depth shown by Figure 13 that it is a non-linear function of time (t) and light intensity. The simple dose-dependence law does not apply to Z_C. Figure 13 also shows that for a given light dose, higher light intensity achieves a lower depth, for large exposure time, with B′t >> 1. This leads to a strategy to use a higher light intensity for faster curing but an extended exposure time (~30%) than what is predicted by BBL is required to achieve the same curing depth as that of lower intensity. This strategy has been clinically demonstrated in corneal cross-link (Lin and Wang, 2016), but not yet in curing of other biomaterials.

The analytic formula of Equation (37) is available only when the oxygen effect is neglected, or K₁₂[O₂] = 0 in Equation (31). The oxygen inhibition effect on Z_C is further investigated by the numerical solutions of Equation (37), including K₁₂[O₂], where Z_C is graphically found by the crossing points of C_EFF = C_T (a threshold value, 0.8) and the curves of C_EFF vs. time for various z. Figure 19 shows that Z_C is a decreasing function of external oxygen supply rate (P′) and the initial oxygen concentration, [O₀], where larger [O₀] or P′ lead to smaller efficacy, thus smaller Z_C due to larger oxygen inhibition effect, as also predicted by Equation (31).

Analysis of Measured Data

Utilization of microfabrication to reduce the deposition steps and to obtain a monolithic product was reported by Alvankarian and Majlis (2015) and Chen et al. (2017), and Wu et al. (2018) in which structures of arrays of pillars in photo-cross-linkable films were measured by irradiation with a periodic array of microscale optical beams under ambient conditions. The optical beams experience a self-focusing non-linearity owing to the photopolymerization-induced changes in refractive index, thereby concentrating light and driving the concurrent, parallel growth of microscale pillars along their path length (Wu et al., 2018). In microfabrication system, formation of radical decreases over depth due to the reduction in light intensity and PS concentration and increase in oxygen inhibition. Under ambient conditions, oxygen diffuses into the film and consumes radicals. The balance of radical production and oxygen inhibition gives rise to the inhibition zone, where the polymerization is completely suppressed.

The curing depth (Z_C) given by the solution of Equation (37) is proportional to the pillar height defined by Chen et al. (2017) where their measured data of their Figure 4 are analyzed as follows. Our Figure 13 shows Z_C vs. time based on Equation (37) that Z_C is proportional to ln (t²I₀), for small t, and reaches its steady-state time (T_S) (when Z_C = 0.9 h), which is scaled by $I_0^{-0.5}$. For example, our calculated T_S = (1,050, 735, 609) s, for I₀ = (5, 10, 15) mW/cm², which shows the similar trend of the measured data of Chen et al. (2017) T_S = (1,050, 600, 480) s. Furthermore, their measured pillar heights ranging 300–600 μm (at t = 600 s), for I₀ = 5–15 mW/cm² are consistent with our Figure 13. Figure 7 of Chen et al. (2017) also showed the similar trend as our scaling law of Z_C vs. light intensity given by Equation (35) that Z_C is scaled by (1/A′)ln(I₀).

Curing depth control is one of the key factors for microfabrication, in which the influencing factors of curing depth include light intensity and exposure time (or dose), the initial concentration of PS and oxygen, and the external supply of oxygen, as shown by Figures 13–19. For a given photo resin thickness of H, Chen et al. (2017) also defined an inhibition zone (Z_N) given by Z_N = H – Z_C. Therefore, a better understanding and manipulation of Z_C (and Z_N) will enable a strong bias differential in the growth between bright (irradiated) and dark regions. The balance of radical production and oxygen inhibition gives rise to the inhibition zone. Figure 8 of Chen et al. (2017) showed a similar trend as our scaling law, in which the log–log plots of inhibition zone vs. light intensity show a linear relationship having slopes depended by various mask spacing ratios.

The above described experimental works have validated certain features of our modeling and the analytic formulas. However, further studies are required to validation other features presented by our modeling, such as (i) the role of PS concentration (C₀) in thick polymers (>1.0 mm), which was assumed as a constant for optically thin polymers (<0.2 mm); (ii) the curing depth at various C₀, as shown by our Figure 16; (iii) dynamic measurement of the light intensity in a strongly depleted PS system, which was assumed as time independent and follows the conventional BLL; however our, theory presented a non-BLL; (iv) as shown by Figure 18, the viscosity effect on the reduction in curing depth; and (v) the role of oxygen external supply rate on the reduction of curing depth, as shown by Figure 19.

FIGURE 18

Figure 18. Curing depth profile. Curing depth vs. time, with (blue) and without (red) the correction of viscosity effect.

FIGURE 19

Figure 19. Curing depth profile. Curing depth vs. time, for various oxygen external supply rates, P′ = (0, 3.0, 5.0, 7.0) × 10⁻⁶ (mM/s), for curves (1, 2, 3, 4); for oxygen initial concentration (A) [O₀] = 0.001 mM and (B) [O₀] = 0.002 mM, for b = 0.0001, C₀ = 0.02 mM and I₀ = 10 mW/cm².

Strategy for Improved Efficacy and Depth

As discussed earlier, a higher light intensity will accelerate the curing speed but suffers a smaller curing depth and efficacy (at steady state). This also indicates that there is a limitation of maximum light intensity and the associate minimum exposure time. Higher PS initial concentration will improve the efficacy, but there is an optimal value (Lin et al., 2019). There are many conventional strategies to reduce oxygen inhibition in photoinduced polymerization: working in an inert or closed environment, increasing the photoinitiator concentration, increasing the light dose or light intensity (for reduced induction time), use of multiple photoinitiators with different rate of initiation, or addition of oxygen scavengers. Chemical mechanisms incorporate additives or suitably functionalized monomers which are insensitive to oxygen, such as the thiol-ene and thiol-acrylate-Michael additive systems (Claudino et al., 2016; Chen et al., 2019a,b). Additive enhancer monomer was proposed to improve the curing (cross-link) efficacy by either reducing the oxygen inhibition effect by stable monomer or increase the lifetime of radicals in clinical applications (Chen et al., 2019b; Wertheimer et al., 2019). The multimonomer system may be applied to industrial materials as an alternative to manipulate the microfabrication. Dual-wavelength (red and UV) photopolymerization was also reported, in which preirradiation of the red light eliminated the oxygen inhibition effect and thus enhanced the conversion efficacy of the UV light (Childress et al., 2019).

Conclusion

Photopolymerization efficacy (C_EFF) is an increasing function of PS concentration (C₀) and the light dose at transient state, but it is a decreasing function of the light intensity, scaled by [C₀/I₀]^0.5 at steady state. The curing depth is an increasing function of C₀ and light dose (time × intensity), but it is a decreasing function of the oxygen concentration, viscosity effect, and oxygen external supply rate. For optically thick polymers, light intensity is an increasing function of time due to PS concentration depletion, which cannot be neglected. Saturation of efficacy profile is governed by the PS depletion rate and light dose. Efficacy is also an increasing function of the effective rate constant, K = k′/k_T^0.5, defined by the radical producing rate (k′) and the bimolecular termination rate. The curing depth is a decreasing function of the oxygen inhibition effect.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

J-TL: concept. D-CC and K-TC: software and data. J-TL and H-WL: supervision. D-CC: financial. J-TL, H-WL, and D-CC: writing and editing. All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest

J-TL was CEO of New Vision Inc.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

D-CC thanks the financial support from the Ministry of Science and Technology (MOST), under grant number MOST108-2221-E-039-003 and China Medical University under grant number CMU 108-MF-98.

References

Alvankarian, J., and Majlis, B. Y. (2015). Exploiting the oxygen inhibitory effect on UV curing in microfabrication: a modified lithography technique. PLoS ONE 10:e0119658. doi: 10.1371/journal.pone.0119658

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J. T., Hudson, S. D., Harrison, C., and Douglas, J. F. (2004). Frontal photopolymerization for microfluidic applications. Langmuir 20, 10020–10029. doi: 10.1021/la049501e

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, F. H., Pathreeker, S. B., and Hosein, I. D. (2017). Synthesis of micropillar arrays via photopolymerization: an in situ study of light-induced formation, growth kinetics, and the influence of oxygen inhibition. Macromolecules 50, 5767–5778. doi: 10.1021/acs.macromol.7b01274

CrossRef Full Text | Google Scholar

Chen, F. M., and Shi, S. (2014). Principles of Tissue Engineering, 4th Edn. New York, NY: Elsevier. doi: 10.1016/B978-0-12-398358-9.00072-0

CrossRef Full Text | Google Scholar

Chen, K. T., Cheng, D. C., Lin, J. T., and Liu, H. W. (2019a). Thiol-ene photopolymerization in thick polymers: kinetics and analytic formulas for the efficacy and crosslink depth. Polymers 11:1640. doi: 10.3390/polym11101640

CrossRef Full Text | Google Scholar

Chen, K. T., Lin, J. T., and Liu, H. W. (2019b). Enhancing radical-mediated photopylomerization efficacy and crosslink depth: kinetic and model of a two-monomer system. Res. Med. Eng. Sci. 8, 853–860. doi: 10.31031/RMES.2019.08.000682

CrossRef Full Text | Google Scholar

Childress, K. K., Kim, K., Glugla, D. J., Musgrave, C. B., and Bowman, C. N. (2019). Independent control of singlet oxygen and radical generation via Irradiation of a two-color photosensitive molecule. Macromolecules 52, 4968–4978. doi: 10.1021/acs.macromol.9b00424

CrossRef Full Text | Google Scholar

Claudino, M., Zhang, X., Alim, M. D., Podgórski, M., and Bowman, C. N. (2016). Mechanistic kinetic modeling of Thiol-Michael addition photopolymerizations via photocaged “superbase” generators: an analytical approach. Macromolecules 49, 8061–8074. doi: 10.1021/acs.macromol.6b01605

PubMed Abstract | CrossRef Full Text | Google Scholar

Cramer, N. B., O'Brien, C. P., and Bowman, C. N. (2008). Mechanisms, polymerization rate scaling, and oxygen inhibition with an ultra-rapid monovinyl urethane acrylate. Polymer 49, 4756–4761. doi: 10.1016/j.polymer.2008.08.051

CrossRef Full Text | Google Scholar

Dendukuri, D., Panda, P., Haghgooie, R., Kim, J. M., Hatton, T. A., and Doyle, P. S. (2008). Modeling of oxygen-inhibited free radical photopolymerization in a PDMS microfluidic device. Macromolecules 41, 8547–8556. doi: 10.1021/ma801219w

CrossRef Full Text | Google Scholar

Fouassier, J. P. (1995). Photoinitiation, Photo-Polymerization, and Photocuring: Fundamentals and Applications; Hanser Gardner Publications. Munich.

Google Scholar

Kotisch, V., Michaudel, Q., and Fors, B. P. (2017). Photocontrolled interconversion of cationic and radical polymerizations. J. Am. Chem. Soc. 139:10665. doi: 10.1021/jacs.7b06661

CrossRef Full Text | Google Scholar

Lin, J. T. (2018). Efficacy S-formula and kinetics of oxygen-mediated (type-II) and non-oxygen-mediated (type-I) corneal cross-linking. Ophthalmol. Res. 8, 1–11. doi: 10.9734/OR/2018/39089

CrossRef Full Text | Google Scholar

Lin, J. T. (2019). Kinetics of enhancer-monomer for corneal cross-linking: two-initiator system. Ophthalmol. Res. 10, 1–6. doi: 10.9734/or/2019/v10i330109

CrossRef Full Text | Google Scholar

Lin, J. T., and Cheng, D. C. (2017). Modeling the efficacy profiles of UV-light activated corneal collagen crosslinking. PLoS ONE 12:e0175002. doi: 10.1371/journal.pone.0175002

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J. T., Liu, H. W., Chen, K. T., and Cheng, D. C. (2019). Modeling the optimal conditions for improved efficacy and crosslink depth of photo-initiated polymerization. Polymers 11:217. doi: 10.3390/polym11020217

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J. T., and Wang, K. C. (2016). Analytic formulas and numerical simulations for the dynamics of thick and non-uniform polymerization by a UV light. J. Polym. Res. 23:53. doi: 10.1007/s10965-016-0934-4

CrossRef Full Text | Google Scholar

O'Brien, A. K., and Bowman, C. N. (2006). Modeling the effect of oxygen on photopolymerization kinetics. Macromol. Theory Simul. 15,176–182. doi: 10.1021/ma051863l

CrossRef Full Text | Google Scholar

Odian, G. (2006). Principles of Polymerization, 4th Edn. New York, NY: John Wiley & Sons, Inc.

Google Scholar

Semchishen, A., Mrochen, A., and Semchishen, V. (2015). Model for optimization of the UV-A/Riboflavin strengthening (cross-linking) of the cornea: percolation threshold. Photochem. Photobiol. 91, 1403–1411. doi: 10.1111/php.12498

PubMed Abstract | CrossRef Full Text | Google Scholar

Wertheimer, C. M., Elhardt, C., Kaminsky, S. M., Pham, L., Pei, Q., Menders, B., et al. (2019). Enhancing rose Bengal photosensitized protein crosslinking in the cornea. Invest. Ophthalmol. Vis. Sci. 60, 1845–1852. doi: 10.1167/iovs.19-26604

PubMed Abstract | CrossRef Full Text | Google Scholar

Wohlers, T., and Caffrey, T. (2016). 3D Printing and Additive Manufacturing State of the Industry. Annual Worldwide Progresss Report, Wohlers Report, Wohlers Associates, Inc.

Google Scholar

Wu, J., Zhao, Z., Hamel, C. M., Mu, X., Kuang, X., Guo, Z., et al. (2018). Evolution of material properties during free radical photopolymerization. J. Mech. Phys. Solid 112, 25–49. doi: 10.1016/j.jmps.2017.11.018

CrossRef Full Text | Google Scholar

Yang, L., Tang, H., and Sun, H. (2018). Progress in photo-responsive polypeptide derived nano-assemblies. Micromachines 9, 296–313. doi: 10.3390/mi9060296

PubMed Abstract | CrossRef Full Text | Google Scholar