Viscoelastic Properties of Polymeric Microneedles Determined by Micromanipulation Measurements and Mathematical Modelling

Microneedles, including dissolvable ones made from biocompatible and biodegradable materials, have been widely studied and can potentially be used for transdermal drug delivery, disease diagnosis (sampling), skin care, etc. Characterizing their mechanical properties is essential, as being mechanically strong enough to pierce the skin barrier is one of the most fundamental and crucial requirements for them. The micromanipulation technique was based on compressing single microparticles between two flat surfaces to obtain force and displacement data simultaneously. Two mathematical models had already been developed to calculate the rupture stress and apparent Young’s modulus, which can identify variations of these parameters in single microneedles within a microneedle patch. In this study, a new model has been developed to determine the viscoelasticity of single microneedles made of hyaluronic acid (HA) with a molecular weight of 300 kDa loaded with lidocaine by using the micromanipulation technique to gather experimental data. The modelling results from the micromanipulation measurements suggest that the microneedles were viscoelastic and their mechanical behaviour was strain-rate dependent, which implies that the penetration efficiency of viscoelastic microneedles can be improved by increasing their piercing speed into the skin.


Introduction
Microneedles are needle-like structures with a length shorter than 1000 µm and a tip width of several tens of microns, and usually, an array of microneedles (e.g., 10 × 10) are fabricated in a patch to be used. They are designed to effectively penetrate the stratum corneum (the topmost layer of skin) without touching the blood vessels or nerves, thus generating no pain and having a low infection risk in contrast with traditional hypodermic needle sticks [1]. They can be used for transdermal drug delivery [2], disease diagnosis (sampling) [3], skin care [4], etc. Microneedles can be classified into several types based on different criteria ranging, from pyramid-shaped, truncated cone-shaped, and other complicated shaped microneedles based on their shapes [2], to glass, metal, and polymer microneedles based on the materials they are made from [5], and to solid, coated, dissolvable (dissolving), and hollow microneedles based on their structure and solubility in body liquid [6]. Dissolvable microneedles made from biocompatible and biodegradable materials that can dissolve or degrade in the human body after administration can potentially have a wide range of applications, including sustained and controlled drug and vaccine delivery and skin care because of their advantages, including no residue of sharp waste after administration and good safety [1,[7][8][9]. For instance, sustained release of a loaded drug is easy to achieve by modifying the dissolvable microneedles with cross-linkable groups [10]. was found to be intrinsic and better to be used to characterise the mechanical strength of microneedles than the rupture force. With the two developed mathematical models, the intrinsic mechanical property parameters as well as their variations for single microneedles can be determined, which can be used to predict their penetration capacity comprehensively and reliably even without knowing their geometrical or structural information. However, it was found that the penetration efficiency of microneedles by impact insertion (high speed) was much higher than by manual insertion (low speed) [23], which implies that the mechanical behaviours of microneedles can vary with the insertion procedure, including the insertion speed, i.e., the microneedles show some speed-dependent behaviour. Therefore, new research needs to be carried out on their viscoelastic and/or viscoplastic behaviour.
In this paper, the elastic model for single microneedles developed in our previous study [8] has been extended to determine their viscoelasticity by accounting for the timedependent behaviours of single microneedles from the experimental data of the micromanipulation tests. HA microneedles loaded with lidocaine have been tested using the micromanipulation technique, and the obtained force versus time/displacement data have been analysed using the developed viscoelastic model to obtain their intrinsic viscoelastic property parameters.

HA Microneedles Loaded with Lidocaine
Microneedles with the shape of a truncated pyramid made from HA and a molecular weight of 300 kDa loaded with lidocaine were used in this study. They were fabricated by micro moulding, and the detailed preparation procedure is reported elsewhere [8]. Briefly, polydimethylsiloxane (PDMS) moulds were first duplicated from a stainless-steel mould and the microneedle body parts were formed by compressing the HA solution (200 mg/mL) dissolved with lidocaine into the PDMS mould using pressured air. After being dried in an anhydrous silica gel environment (room temperature and humidity) for 0.5 h, the microneedle base plate was prepared by adding pure hyaluronic acid (HA) to the PDMS mould. Finally, the microneedle patches were peeled off from the PDMS mould after drying in an anhydrous silica gel environment (room temperature and humidity) for another 4 h. Their formula is illustrated in Table 1. Table 1. HA microneedles loading with lidocaine.

Micromanipulation of Single Microneedles
The micromanipulation technique is based on the diametrical compression of single microparticles/microneedles between two parallel flat surfaces to different deformations, including rupture, during which force-displacement/time data are obtained simultaneously. The detailed principle has been reported elsewhere [24][25][26][27]. In the micromanipulation tests, the bottom of a microneedle patch was fixed to the surface of a glass slide using doublesided tape (30401, Deli Group Co., Ltd., Ningbo, China). The glass slide was fixed onto the sample stage of the micromanipulation rig. Single microneedles of Lido-HA300kDa were compressed to a certain deformation and held for 10 s under the transducer probe, as illustrated in Figure 1. The ambient temperature and humidity were 17 ± 2 • C and 38 ± 3% Rh, respectively. transducer probe, as illustrated in Figure 1. The ambient temperature and humidity were 17 ± 2 °C and 38 ± 3% Rh, respectively. In order to assess whether the double-sided tape can generate any viscoelastic behaviour in the system, a small glass with a similar dimension (6.5 mm × 6.5 mm) to the microneedles' base was placed on top of the double-sided tape, which had a thickness of ~ 0.1 mm (see Figure S1 in the supplementary information). The diameter of the probe was around 60 μm. The compression displacement of the force transducer probe on the small glass slide was set to be 8-30 μm, and the corresponding force versus sampling time was recorded.

Derivation of the Viscoelastic Model
The schematic diagram of a single microneedle under compression is illustrated in Figure 2. Assuming the microneedle is perfectly sharp initially, the relationship between the force F and displacement δ of single microneedles under compression within the elastic limit can be expressed using the elastic model in Equation (1) developed in the previous study [8].
where is Young's modulus of the microneedle. From similarity, ℎ′ ℎ ℎ′ (2) can be calculated using the following equation for the microneedles studied in this research with a shape of a truncated pyramid:

4
(3) The mathematical relationship between Young's modulus and the shear modulus can be expressed using Equation (5) [28].  In order to assess whether the double-sided tape can generate any viscoelastic behaviour in the system, a small glass with a similar dimension (6.5 mm × 6.5 mm) to the microneedles' base was placed on top of the double-sided tape, which had a thickness of~0.1 mm (see Figure S1 in the supplementary information). The diameter of the probe was around 60 µm. The compression displacement of the force transducer probe on the small glass slide was set to be 8-30 µm, and the corresponding force versus sampling time was recorded.

Derivation of the Viscoelastic Model
The schematic diagram of a single microneedle under compression is illustrated in Figure 2. Assuming the microneedle is perfectly sharp initially, the relationship between the force F and displacement δ of single microneedles under compression within the elastic limit can be expressed using the elastic model in Equation (1) developed in the previous study [8].
where E is Young's modulus of the microneedle. From similarity, A 1 can be calculated using the following equation for the microneedles studied in this research with a shape of a truncated pyramid: Combination of Equations (1)-(3) leads to The mathematical relationship between Young's modulus E and the shear modulus G can be expressed using Equation (5) [28].
where ν is the Poisson's ratio. Letting and combining Equations (4)-(6) leads to the following equation: where is the Poisson's ratio. Letting 4 ℎ (6) and combining Equations (4)-(6) leads to the following equation: Figure 2. Schematic diagram of a single microneedle under compression. is the area of the microneedle bottom, ℎ is the initial length of the microneedle body, and are the initial half side length of a quadrangular microneedle base and tip respectively, and ℎ is the height of the missing tip.
Following the approach of [28], replacing the constant term 2 with the time-varying shear modulus leads to Viscoelasticity can be seen as a combination of a viscos unit and a linear elastic unit, which are most commonly described by a Newtonian dashpot and a Hookean spring, respectively [29]. Several models can be used to describe the viscoelastic behaviour of materials, and the Maxwell and Kelvin models are the two most basic ones, which can be seen as the serial and parallel combination of the viscos and linear elastic units. Other complex viscoelastic models can be built from different combinations of Maxwell and Kelvin models [26,29]. The generalised Maxwell model (also called the Prony series [28]), which is built by the combination of n terms of the basic Maxwell model, is more general and semi-empirical [26]. It has been successfully used to describe the spherical indentation and relaxation of soft biological tissues [30] and the compression and relaxation of agarose micro-particles [28]. Thus, it is also used in this paper to model the relaxation of microneedles, which are made from biocompatible materials. During the relaxation, the time-dependent shear modulus can be expressed using the following equation [28]: (9) ℎ ℎ initial state compression at Following the approach of [28], replacing the constant term 2G with the time-varying shear modulus G(t) leads to Viscoelasticity can be seen as a combination of a viscos unit and a linear elastic unit, which are most commonly described by a Newtonian dashpot and a Hookean spring, respectively [29]. Several models can be used to describe the viscoelastic behaviour of materials, and the Maxwell and Kelvin models are the two most basic ones, which can be seen as the serial and parallel combination of the viscos and linear elastic units. Other complex viscoelastic models can be built from different combinations of Maxwell and Kelvin models [26,29]. The generalised Maxwell model (also called the Prony series [28]), which is built by the combination of n terms of the basic Maxwell model, is more general and semi-empirical [26]. It has been successfully used to describe the spherical indentation and relaxation of soft biological tissues [30] and the compression and relaxation of agarose microparticles [28]. Thus, it is also used in this paper to model the relaxation of microneedles, which are made from biocompatible materials. During the relaxation, the time-dependent shear modulus can be expressed using the following equation [28]: Similarly, the corresponding force can be expressed as [28]. where τ i , i = 1, 2, . . . , n is the relaxation time. The relationship between the coefficients C 0 , C i and B 0 , B i are and where δ 0 is the maximum displacement corresponding to the end of loading, and Φ i is the ramp correction factor given by [28] where t R is the rising time from the onset of loading to the onset of relaxation.
It is worth pointing out that the value of the relaxation times τ 1 , τ 2 , · · · indicate how fast each relaxation term takes place and the magnitudes of B 1 , B 2 , · · · suggest the significance of each relaxation term. Smaller values of relaxation times mean faster relaxation, and bigger values of magnitudes have a more significant influence on overall relaxation.
The instantaneous shear modulus G 0 and long-term shear modulus G ∞ can then be estimated from the coefficients C 0 , C 1 , C 2 , . . . , C n [28]. and The corresponding instantaneous Young's modulus, E 0 , and long-term Young's modulus, E ∞ , can be calculated by [28]. and

Determination of Viscoelastic Parameters Based on Numerical Optimization
The viscoelastic model with degree n ≥ 2 as shown in Equation (10) cannot be used analytically to determine the model parameters from experimental data, but the viscoelastic parameters can be obtained by numerical optimization to minimize the sum of the squared difference between the predicted forceF i and the experimental force F i based on the following equation: where the predicted force,F i , is calculated using Equation (10) and N is the number of data points during the holding experiment.

Prediction of the Loading Force Using the Viscoelastic Parameters
With the parameters obtained from the force relaxation data [30], the following equation can be used to predict the loading force data: If the compliance of the force transducer is negligible, δ(t) can be expressed as [28] δ(t) = Vt (20) where V is the moving speed of the micromanipulator holding the force transducer. Substituting Equations (9) and (20) into Equation (19) leads to the following equation: which can be used to predict the force during the compression of single microneedles if the compliance of the force transducer is negligible. However, if the compliance of the force transducer is not negligible relative to the imposed displacements, significant errors can result from using Equation (21). Here, two strategies can be used. One is using the average moving speed of the force transducer tip, which can be determined from the corresponding displacement versus time data, to substitute the moving speed of the micromanipulator V. Another is using numerical solutions. The displacement of the force transducer tip δ(t) is given by [27] where c is the compliance of the force transducer. Substituting Equation (22) into Equation (19) results in The second part, which includes F(t) itself, cannot be solved mathematically. To deal with this problem, a numerical solution can be used by decomposing the loading history into a sum of infinitesimal loading steps based on the Boltzmann superposition principle [31] where, ∆t i , i = 1, 2, · · · , is the infinitesimal time interval in which step loading is applied. In practice, the acquisition time (t s ) can be used as ∆t i , i = 1, 2, · · · to decompose the loading history, then Equation (24) can be written as follows: where j is the sampling number. Equation (25) can also be written as and from which the following equation can be derived The corresponding discrete form of Equation (22) can be written as and Thus A combination of Equation (28) and Equation (31) leads to the following recursive equation and Equations (32) and (33) can be used to predict the forces during loading in the typical loading-holding experiment of microneedle tests in a numerical recursive way without the restriction of negligible force transducer compliance as required when Equation (21) is used.

Force-Time Data of Lido-HA300kDa Microneedles under Compression and Holding
Typical curves of force versus time data during compression and holding of the glass slide (as control) and single Lido-HA300kDa microneedles whose bases were placed on the double-sided tape are illustrated in Figures 3a and 3b, respectively. Clearly, the force increased during compression in both cases. During the holding of the glass slide (which is considered to be very rigid) on the tape, there was no significant force relaxation, whereas the force dropped notably during the holding of the microneedle, which suggests that some viscoelastic behaviour of the microneedle existed.

Elastic Analysis
If the viscos behaviour is neglected, the force-displacement data of the Lido-HA300kDa microneedle during compression can be fitted into the elastic model in Equation (4) as illustrated in Figure 4. The coefficient of determination ( ) is 0.99, which suggests good fitting. The Youngs modulus obtained is 422.8 MPa.

Elastic Analysis
If the viscos behaviour is neglected, the force-displacement data of the Lido-HA300kDa microneedle during compression can be fitted into the elastic model in Equation (4) as illustrated in Figure 4. The coefficient of determination (CoD) is 0.99, which suggests good fitting. The Youngs modulus obtained is 422.8 MPa.
(a) (b) Figure 3. Typical force-time curves of the loading-holding test of the glass slide (a) and Lido-HA300kDa microneedles (b) (tip width 2r = 16.1 μm, bottom width 2r = 250 μm, h = 700 μm) whose bases were placed on the double-sided tape. A: Onset of compression. B: Onset of holding (end of compression). C: End of holding. Compression speed was 2 μm/s.

Elastic Analysis
If the viscos behaviour is neglected, the force-displacement data of the Lido-HA300kDa microneedle during compression can be fitted into the elastic model in Equation (4) as illustrated in Figure 4. The coefficient of determination ( ) is 0.99, which suggests good fitting. The Youngs modulus obtained is 422.8 MPa.

Viscoelastic Analysis
Letting the degree = 2 and fitting the force-time data during the holding of the microneedle in Figure 3b into Equation (10) optimised with Equation (18) using Microsoft ® Excel Solver, the coefficients , , and the time constants and can be obtained as shown in Table 2. It can be seen that the viscoelastic model with two relaxation times can predict the experimental holding data of the microneedle in Figure 3b Table 3 obtained using Equations (11) and (12).

Viscoelastic Analysis
Letting the degree n = 2 and fitting the force-time data during the holding of the microneedle in Figure 3b into Equation (10) optimised with Equation (18) using Microsoft ® Excel Solver, the coefficients B 0 , B 1 , B 2 and the time constants τ 1 and τ 2 can be obtained as shown in Table 2. It can be seen that the viscoelastic model with two relaxation times can predict the experimental holding data of the microneedle in Figure 3b well (CoD = 0.98), as illustrated in Figure 5. The rising time t R = 2.60 s, the displacement up to the end of compression δ 0 = 2.52 µm, and the ramp correction factors Φ 1 and Φ 2 are 1.43 and 384.88, respectively. The coefficients C 0 , C 1 , C 2 are shown in Table 3 obtained using Equations (11) and (12). Table 2. Calculated values of the coefficients and relaxation times of the microneedle in Figure 3b.

Width (µm)
B 0 (mN) B 1 (mN) B 2 (mN) τ 1 (s) τ 2 (s) 16  The relaxation time of the second relaxation (0.32 s) is nearly one order smaller than the loading time (2.60 s), which indicates the relaxation is quite rapid. Although the first relaxation time (3.82 s) is in the same order as the loading time (2.60 s), the magnitude of the first relaxation ( = 1.55 mN) is much smaller than the second one ( = 534.35 mN), which suggests that the second relaxation is dominant for the relaxation of the tested microneedle in Figure 3b.  The relaxation time of the second relaxation (0.32 s) is nearly one order smaller than the loading time (2.60 s), which indicates the relaxation is quite rapid. Although the first relaxation time (3.82 s) is in the same order as the loading time (2.60 s), the magnitude of the first relaxation (B 1 = 1.55 mN) is much smaller than the second one (B 2 = 534.35 mN), which suggests that the second relaxation is dominant for the relaxation of the tested microneedle in Figure 3b.
The instantaneous shear modulus G 0 , long-term shear modulus G ∞ , instantaneous Young's modulus E 0 and the long-term Young's modulus E ∞ of the microneedle in Figure 3b have been obtained using Equations (14)- (17), and the values are shown in Table 4. The Poisson's ratio ν of 0.5 was used. In last section, the Youngs modulus value without considering the viscos behavior is 422.8 MPa, which is expected to be between the values of the instantaneous Young's modulus (513.3 MPa) and long-term Young's modulus (343.2 MPa). Other tested Lido-HA300kDa microneedles yielded comparable results.
The force data during loading for the microneedle in Figure 3b have been predicted using the parameters obtained from the holding (i.e., force relaxation) data, which are shown in Figure 6. It can be seen that the results obtained from the approach with the average moving speed of the force transducer tip (F-P2) and the numerical approach (FP-1) are both consistent with the experimental data, and thus they both can be used to predict the compression force utilising the viscoelastic parameters obtained from the data of the force relaxation corresponding to the holding. The predicted force F-P3 deviates from the experimental data significantly as the compliance of the transducer c is 0.5 µm/mN, which is not negligible with such a big force in the order of mN or above.
Materials 2023, 16, x FOR PEER REVIEW 11 of 14 Figure 6. Predicted force data in comparison with the experimental data corresponding to compression of the microneedle. F is the experimental force. F-P1 is predicted by the numerical method using Equation (32). F-P2 is predicted using Equation (21) with the average moving speed of the force transducer tip (= ⁄ ). F-P3 is obtained using Equation (21) with the pre-set moving speed of the micromanipulator holding the force transducer.
Totally, twenty-six single Lido-HA300kDa microneedles were tested. Their experimental data during compression were analysed using the elastic model by neglecting the viscos behaviour and the obtained mean Young's modulus value was 496.5 ± 31.6 MPa, which lies between the obtained mean instantaneous and long-term Young's modulus values (619.0 ± 37.2 MPa and 390.5 ± 24.1 MPa, respectively) obtained from the analysis of their holding data using the developed viscoelastic model. The variations of the viscoelastic property parameters, including the instantaneous and long-term Young's moduli,  Figure 6. Predicted force data in comparison with the experimental data corresponding to compression of the microneedle. F is the experimental force. F-P1 is predicted by the numerical method using Equation (32). F-P2 is predicted using Equation (21) with the average moving speed of the force transducer tip (=δ 0 /t R ). F-P3 is obtained using Equation (21) with the pre-set moving speed of the micromanipulator holding the force transducer.
Totally, twenty-six single Lido-HA300kDa microneedles were tested. Their experimental data during compression were analysed using the elastic model by neglecting the viscos behaviour and the obtained mean Young's modulus value was 496.5 ± 31.6 MPa, which lies between the obtained mean instantaneous and long-term Young's modulus values (619.0 ± 37.2 MPa and 390.5 ± 24.1 MPa, respectively) obtained from the analysis of their holding data using the developed viscoelastic model. The variations of the viscoelastic property parameters, including the instantaneous and long-term Young's moduli, can mainly result from the physiochemical properties and composition of the polymer matrix used to prepare the microneedles, which has been verified in [8], where variations of the mechanical strength over a patch of the four tested microneedle samples made from HA with and without loaded drugs were identified. It was reported that the value of Young's modulus of human skin also depends on the loading speed and direction as well as the test methods, which were found to be 5-100 kPa by indentation tests, 0.025-140 MPa by tensile and torsion tests, and 25-260 kPa by suction tests [32]. The instantaneous and long-term Young's modulus values of Lido-HA300kDa microneedles are significantly greater than those of human skin, and therefore they could be strong enough to reliably pierce it.
Moreover, the difference between the instantaneous and long-term moduli is quite significant, implying that different penetration efficiency can be obtained by adjusting the insertion speed, i.e., with a higher insertion speed, more efficient penetration can be achieved because the apparent Young's modulus value should increase with compression speed.

Further Discussion
To precisely determine the intrinsic mechanical property of single microneedles, the micromanipulation technique had been used, and an elastic model has been developed in our previous studies [8] to determine Young's modulus of single microneedles from the experimental data obtained at a given compression speed, which may be considered as the apparent Young's modulus. Young's modulus is an intrinsic mechanical property of microneedles that are independent of their geometrical or structural information. However, the microneedles made from dissolvable polymers usually show some speed-dependent (viscoelastic) behaviour, and their mechanical properties, including the apparent Young's modulus, can vary with the strain rate used to test them. Therefore, in this paper, the elastic model has been extended to analyse the viscoelastic behaviours of microneedles, and a mathematical model has been developed to determine their viscoelastic parameters, including the instantaneous and long-term Young's moduli. The apparent Young's modulus value obtained from using the purely elastic model [8] is a value at a certain compression speed (i.e., 2 µm/s), whereas the instantaneous and long-term Young's moduli obtained from using the viscoelastic model are the upper and lower limits of Young's modulus values, which correspond to compressions at infinitely large and infinitely small speeds, respectively. Results of Lido+HA300kDa microneedles show that their Young's modulus values at a compression speed of 2 µm/s (obtained from the elastic model) were between their instantaneous and long-term Young's modulus values, and the difference between the latter two is quite significant, which suggests that the microneedles show higher stiffness when they penetrate the skin with an increasing insertion speed, and consequently a higher penetration efficiency can be obtained. Although it was found that drug loading could decrease the mechanical strength of microneedles [8,10], reliable piercing can still be achieved by increasing the insertion speed using a mechanical insertion device, such as an impact-insertion applicator [23] for microneedles loaded with drugs.

Conclusions
In conclusion, a viscoelastic model has been developed in this paper to determine the viscoelastic parameters of single microneedles tested using the micromanipulation technique. In combination with the rupture strength model and elastic model developed in the previous study [8], the intrinsic mechanical property parameters as well as their variations in single microneedles can be precisely determined from the micromanipulation tests, which can be used to predict their penetration capacity comprehensively and reliably even without knowing their geometrical or structural information. A numerical recursive equation (Equation (32)) has been developed to predict the loading force using the viscoelastic parameters obtained from the holding data, which has the advantage that the compliance of the force transducer is included in comparison with the analytical method in which the compliance was negligible [28].
Moreover, the three mathematical models can also be used to analyse the experimental data of microneedle arrays tested using other techniques, including the displacement-force test station [15] and the micro-mechanical test machine with a microforce sensor [18]. Therefore, they can find wide applications in academic research and the development of new products based on microneedles.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ma16051769/s1, Figure S1: Loading-holding test of the glass slide whose base was placed on the double-sided tape.

Conflicts of Interest:
The authors declare no conflict of interest.