A design of experiments approach to identify the influencing parameters that determine poly-D,L-lactic acid (PDLLA) electrospun scaffold morphologies

Electrospun fibrous materials have increasing applications in regenerative medicine due to the similarity of fibre constructs to the morphology of certain extracellular matrices. Although experimentally the electrospinning method is relatively simple, at the theoretical level the interactions between process parameters and their influence on the fibre morphology is not yet fully understood. Here, we hypothesised that a design of experiments (DoE) model could determine combinations of process parameters that result in significant effects on poly-D,L-lactic acid (PDLLA) fibre morphology. The process parameters used in this study were applied voltage, needle-to-collector distance, flow rate and polymer concentration. Data obtained for mean fibre diameter, standard deviation (SD) of the fibre diameter (measure of fibre morphology) and presence of ‘beading’ on the fibres (beads per μm2) were evaluated as a measure of PDLLA fibre morphology. Uniform fibres occurred at SDs of ≤500 nm, ‘beads-on-string’ morphologies were apparent between ±500 and 1300 nm and large beads were observed at ±1300–1800 nm respectively. Mean fibre diameter was significantly influenced by the applied voltage and interaction between flow rate and polymer concentration. Fibre morphology was mainly influenced by the polymer concentration, while bead distribution was significantly influenced by the polymer concentration as well as the flow rate. The resultant DoE model regression equations were tested and considered suitable for the prediction of parameters combinations needed for desired PDLLA fibre diameter and additionally provided information regarding the expected fibre morphology.


Introduction
Electrospinning is widely used in the manufacture of biologically relevant scaffolds for cell culture, enabling the formation of micro or nanoscale fibres and porous matrix structures [1]. Scaffolds produced by this method exhibit similarities in physical and material properties to the native extracellular matrix (ECM), facilitating their application in regenerative medicine and tissue engineering as in vitro culture platforms for mammalian cells [1][2][3][4][5][6][7][8][9]. In the electrospinning process, a polymer solution or melt is accelerated through a charged needle towards a grounded or opposite polarity collector. A high voltage is applied to the needle, which acts on the meniscus solution droplet at the tip. The rounded droplet extends to form a conical shape known as a Taylor cone [10], in which the electric field exceeds the surface tension of the droplet resulting in a fibre jet being ejected from the Taylor cone towards the collector. Solvent evaporation occurs while the fibre jet travels through the atmosphere, leading to solid polymer fibre deposition on the collector [11]. Flat or rotating mandrel collectors are commonly used [9]. Although electrospinning can produce well-defined fibres, process parameters can greatly affect fibre morphology, and these factors are poorly understood. Many parameters may change the resultant fibre morphologies; these factors include applied voltage, collector plate distance from the needle, polymer solution concentration, solvent type and flow rate. The interaction between these process parameters is complex, with variation of one factor often altering another, as shown in figure 1 [12][13][14][15][16][17][18][19][20]. It is therefore challenging to investigate each individual parameter experimentally and the resultant procedure is time consuming. For reproducible electrospinning to occur, the key process parameters must be identified and stably controlled.
In general, three fibre morphologies can arise from electrospinning; uniform fibres, 'beads-on-strings' and 'big beads' (figure 2). The occurrence of beading is thought to be as a result of flow rate or polymer concentration [21][22][23], with beading impacting upon the reproducibility, biological significance and downstream cell viability [24]. As a result, much effort has been invested to identify process parameters that will yield uniform fibres. Although bead formation is a common phenomenon, the understanding of the parameter combinations resulting in beaded morphologies is still unclear. Electrospinning process optimisation is a common topic in published literature as researchers attempt to understand the specific effect of multiple variables [19,21,[25][26][27][28][29].
Much research has been reported on the impact of such parameters (on fibre diameter for example), however little is known about the interaction between these parameters (i.e. how they influence each other). To date, experimental design towards understanding process parameter effects has mostly involved one factorial approaches (see more info in the supplementary data; figure S1 is available online at stacks.iop.org/ BMM/12/055009/mmedia). This singular focus excludes the interrelationship between interacting parameters, and as mentioned above, results in time consuming optimisation experiments. Design of experiments (DoE) or quality by design (QbD) approaches use fractional design to identify the most influential parameters. This method takes into account the interaction of parameter interrelationship, using linear regression and analysis of variance (ANOVA) mathematical models [30], identifying the significantly impacting factors and parameter interactions on the resultant fibre characteristics. Different combinations of the factors (process parameters such as applied voltage etc) levels (high or low values of the parameter settings, for example applied voltage 10, 15 and 20 kV) are tested and resultant outputs (such as fibre diameter etc) are used to define factor combination interactions. Building on the understanding of parameter interaction and interdependence, the DoE approach contributes to the overall understanding of the electrospinning process as a manufacturing technique. Currently there is limited research focusing on DoE approaches in electrospinning [31][32][33][34][35][36][37]. Desai and Sung [38] described a DoE approach using a two level fractional design of four factors. Polymer concentration and needle-collector distance, individually and in combination, were observed to have a significant impact on the fibre thickness and bead density (number of beads) of fibres fabricated from polyaniline/ poly (methyl methacrylate) (PANI/PMMA). Cui et al [39] used a three level fractional design of six factors to investigate the effect of parameters (voltage, polymer conc., molecular weight (MW), solvent system, flow velocity, and nozzle size) on the mean fibre diameter and percentage of beading of PDLLA fibres. A significant impact of MW and concentration of the polymer on the fibre diameter was observed. Likewise, a significant effect of MW, concentration, and solvent system used, on the percentage of beading was identified. A more recent and larger DoE investigation by Seyedmahmoud et al [40] used different polymer MWs as a variable while keeping the polymer concentration constant (12% w/v), to fabricate bead-free poly-L-lactic acid (PLLA) fibres. This DoE was proposed by the authors to be more robust due to the size of the model: five factors, three levels and replicated runs (96). A significant impact of polymer MW on the fibre diameter was observed. To date, there are few published studies using DoE models to identify parameters that would yield uniform fibre formation with little to no beading. When embarking upon a new DoE approach for electrospinning, an initial optimisation phase is performed to identify the key parameter values to be selected for the model. Most chosen parameters are selected based on the production of uniform fibres whilst excluding the 'undesirable' morphologies of bead formation that can be observed in the electrospinning process. In adopting this 'trial and error' approach, there is a lack in understanding of the impact of process parameters on the different fibre based morphologies. Investigating these features (such as 'big beads' and 'beads-on-strings'; figures 2(A)/ (B)), in addition to uniform fibres, would provide a better understanding of which parameters lead to which morphologies. This should reduce optimisation time and give an indication of the most appropriate process parameter window required to achieve the desired fibre morphology. A DoE model can also provide a linear regression equation that can be used to predict the fibre diameter and the chance of bead formation within the chosen set of parameters.
In this paper therefore, we detail a DoE method focused on finding the parameters which significantly influence the fibre diameter, morphologies and bead distribution observed in electrospinning of PDLLA fibres. A simple three level design of four factors and three input levels (high; 1, middle; 0 and low; −1) in combination was performed. This resulted in a large dataset, which has the benefit of increased accuracy of statistical estimates and confidence intervals (CIs), which in turn are crucial for a systematic and comprehensive investigation of electrospinning parameters [40].   70:30) was purchased from Evonik Industries, Germany, Essen.

DOE model
A DOE approach of four factors and three levels was used. From 25 different process parameters combinations, 81 scaffolds, including triplicate samples and nine replicates of middle point, were generated. A 'randomisation strategy' was used to reduce variability. In addition, a one single operator approach was used to avoid un-necessary variability in the model. Linear regression, ANOVA and data analysis were performed by MODDE 10.1.1 and Minitab 17.2.1. The following linear regression equation (equation (1)) was used to identify the most significant impacting parameters (X 1 , X 2 , X 3 , X 4 ) on the fibre diameter, standard deviation (SD) and beads per μm ( ) The linear regression output was fitted for each Y output of the experimental data to obtain an 'ordinary least square' linear estimation for Y(X 1 , X 2 , X 3 , X 4 ). This is a third order complete model with X i , X j , X k , and X l the independent electrospinning parameters (applied voltage, distance, flow rate and polymer concentration) that influence the response Y i (mean fibre diameter, SD or number of beads per μm 2 ). β 0 is the model constant; β i is the linear coefficient representing the influence factor of individual independent parameters on Y i . β ij is the cross first order interactions, product coefficients representing the influence factor of two independent parameter interactions on Y i . β ijk is the second order interactions, product coefficients representing the influencing factor of three independent parameter interactions on Y i . β ijkl is the cross third order interaction coefficient representing the influencing factor of all independent parameter interactions on Y i and ε is the error. The higher the β coefficient, the more influence it has on Y i . Using this approach, the most important individual or combinations of independent parameters can be identified.

Electrospinning
Electrospun solutions were prepared, the day before the planned experiment, in Ac:DMF (1:1) with the desired concentration of PDLLA at ambient temperature by stirring. Obtained scaffolds were fabricated using a Climate Control EC-CLI platform with rotating target collector EM-RTC electrospin unit (IME Technologies, the Netherlands, Geldrop) with a 0.8 mm internal diameter (21 gauge) flat tip needle. The conducting surface was reduced to a width of 5 cm (approx. 100 cm 2 surface area) by covering either side of the rotatating mandrel (7.5 cm diameter) with parafilm. The climate control system was started 30 min before the run to obtain constant temperature and humidity of 20°C and 44% respectively. A total volume of 2 ml polymer solution per run was used in the spinning process for each scaffold batch. The tubing and needle were carefully rinsed with acetone to remove any of the previous samples prior to every electrospinning run.

Scanning electron microscopy (SEM)
For each obtained scaffold a 1×1 cm sample was cut out of the centre, mounted on a pin stubs (12.5 mm) with carbon tabs (12 mm) and gold sputter coated for 300 s at 25 mA (EM SCD005, Leica Microsystems). All scaffolds were analysed by SEM (JSM-6060LV, Jeol LtD Japan) at ×1000, ×2000, and ×5000 magnifications (mag.) with six different fields of view captured for each scaffold.

Fibre diameter, SD and number of beads per μm 2 measurements
Fibre diameter measurements were performed with ImageJ 1.48v software on the six ×5000 mag. images; 20 measurements per ×5000 mag. SEM images (120 in total). Beads were included in the measurements and measured on the widest point of the bead. Mean fibre diameter and SD of the measurements were calculated by GraphPad Prism software. Number of beads per μm 2 measurements were performed with ImageJ 1.48v software on three ×2000 magnification SEM images. Beads per image were counted and beads per μm 2 were calculated.

Results
In this paper a DoE approach was performed to identify the most significant electrospinning parameter or parameter interaction effects on the formation of fibres, their morphology and amount of beads per μm 2 . A preliminary study, supplementary table S1, identified the sets of parameter levels (table 2) for the DoE model used to produce fibres of all three morphologies. From these sets, a DoE model of four factors, three levels and three replications was considered suitable, resulting in 81 electrospinning runs in total (25 different parameter sets, three replicates for each parameter set, nine replicates of the middle point).

Fibre diameter (Y 1 )
A normal distribution was observed throughout the data with a mean of 871±189 nm ( figure 3(A)). Linear regression was observed (the corresponding linear regression equation can be found in supplementary table S3, equation (A.1)). A pareto plot ( figure 3(B)) showed applied voltage as the most significant influencing parameter (X 1 ), followed by the flow rate and polymer concentration combined (X 3 X 4 ) on the mean fibre diameter and the polymer concentration on its own (X 4 ). The corresponding linear regression equation can be simplified to include only the significant impacting factors in the regression equation, as the non-significant impacting factors will generate a small change in the outcome Y of the equation due to their low β factor (supplementary table S3, equation (S2)). The 95% CI for estimates was calculated giving a 95% CI of ±32 nm for a full prediction and a 95% CI of ±35 nm for the simplified prediction.

Fibre morphology (Y 2 )
SDs in the range of ±200 to ±1800 nm were observed throughout the DoE. Figure 3(C) shows ±400-500 nm and ±1000-1200 nm as the most frequently observed SDs. SEM images with SD ranges 500 nm showed scaffolds with uniform fibres.
Beads-on-strings formation was observed between 500 and 1300 nm and large beads were observed within the range of ±1300-1800 nm. Therefore, SD was concluded to be an appropriate measure of fibre morphology. A pareto plot showed a significant influence of only the polymer concentration (X 4 ; figure 3(D)) on the corresponding fibre morphology. Linear regression was observed (the corresponding full and simplified linear regression equations can be found in supplementary table S3, equations (S3) and (S4)); 95% CI was calculated giving a 95% CI of 359 nm for a full prediction and a 95% CI of 358 nm for the simplified prediction.

Beads per
The number of beads per μm 2 had a direct correlation with the SD ranges observed for the fibre diameters. More beads per μm 2 were observed at higher SDs for big bead morphologies (±1-3 beads per μm 2 ), the beads-on-string morphology correlated to <1 beads per μm 2 , while uniform fibres correlated with 0-0.2 beads per μm 2 (figure 3(E)). No beading was observed at the highest polymer concentration 25% w/v PDLLA at the lowest flow rate (1 ml h −1 ). The highest number of beads (3 beads per μm 2 ) was observed for the lowest concentration 15% w/v PDLLA and the highest flow rate of 5 ml h −1 . The pareto plot showed a significant impact of polymer concentration (X 4 ) and flow rate (X 3 ) on the number of beads per μm 2 ( figure 3(F)).

Test model predictions
The predictive equations (supplementary table S3, equations (S1)-(S6)) were tested by choosing different random combinations of parameters within the boundary of the DoE model (table 1) and testing the resulting prediction against measurements generated from scaffolds fabricated using those test parameters. Mean fibre diameters from the obtained scaffolds were within the 95% CI of the predicted mean fibre diameters for all test runs (figure 4(A)). A much higher 95% CI was observed for the SD and number of beads per μm 2 . The SD of run B/C and the number of beads per μm 2 of all runs were within the 95% CI while the SD of run A was further away from the model (figures 4(B)/(C)).

Discussion
We aim to show that a DoE approach can identify the most significantly impacting parameters on the resulting PDLLA fibre morphology. As the difference in isomeric properties of PLA has been investigated in the literature [41][42][43], we have exclusively focused this work on PDLLA polymer. The preliminary study identified a correlation between the variation of fibre diameters over the different morphologies observed in the electrospinning process. These three different morphologies were termed 'big beads', 'beads-onstrings' and 'uniform fibres' (figure 2). Although the SD gives an indication of morphology, it does not give a measure of the number of beads. Therefore, number of beads per μm 2 was included in the DoE model. The significant parameters that influenced mean fibre diameter, SD and number of beads per μm 2 were identified by the DoE model. Parameter combinations and corresponding average fibre diameter, SDs and    ; table S2. For a stable electrospinning process, the force (flow rate) required to push the polymer solution into the Taylor cone must be equal to the force generated by the jet pull rate that carries the solution away. This jet pull rate is created by the difference in electrical potential [44]. A good balance between the flow rate and the applied voltage is therefore important [45] ( figure 5(A)). Polymer concentration has a direct impact on the whipping and stretching of the polymer jet. A decrease in whipping occurs with increasing polymer concentration due to the physical weight of polymer and polymer/solvent interaction within the jet. This results in an increase in fibre diameter and explains the significant interacting effect of flow rate and polymer concentration on the resultant fibre diameter observed here as well as the significant impact of applied voltage ( figure 3(B)).
The three observed morphologies corresponded to particular SD ranges due to including the bead diameters in the fibre diameter measurements. By including the beads in the measurement the SD increased when the scaffold was heterogeneous in diameter (due to this beading) and decreased when the fibres were uniform and more homogenous. Therefore, lower SDs corresponded to uniform fibres (homogenous morphologies) while larger SDs corresponded to large bead formation (heterogeneous morphologies). An increase in beading was observed when the PDLLA concentration was reduced to both 20 and 15% w/v. The DoE results confirmed this observation showing that the fibre morphology was only significantly affected by the polymer concentration (X 4 ; figure 3(D)). The number of beads per μm 2 value had a direct correlation with the SD ranges observed for the fibre diameters. The amount of beading was observed to increase at the lowest concentration of PDLLA and the highest flow rate. Whilst no beading was observed for higher concentrations of PDLLA and lower flow rates. This observation was confirmed by the DoE output showing a significant impact of polymer concentration (X 4 ) and flow rate (X 3 ) on the number of beads per μm 2 ( figure 3(F)). This we attributed to the instability of the Taylor cone, which increased when a high flow rate was used with a low polymer concentration solution [44] ( figure 5(B)). The number of beads per μm 2 was therefore directly correlated with the flow rate used.
The obtained fibre diameters of the test runs were within the 95% CI of the predicted mean fibre diameters for all test runs ( figure 4(A)). A much higher 95% CI was observed for the SD and number of beads per μm 2 . This was due to the inability to control the unstable Taylor cone, which made the prediction of the exact size of the beads and number of beads per μm 2 less accurate. A beads-on-strings formation was observed (figures 4(D)-(F)) when morphology SD predictions fell within the beads-on-strings SD ranges (table 3). The number of beads per μm 2 data showed the same correlation of the predicted number of beads ranges to SD ranges. We considered a more appropriate way to predict fibre morphology and beading was to group the calculated SD or number beads per μm 2 to the corresponding morphology, i.e. the extent of beading (table 3).
The concentration of the polymer solution significantly affects Taylor cone stability, as observed in the DoE, with low polymer concentrations leading to jet instability. This is because polymer concentration influences the viscosity of the solution which in turn affects Taylor cone stability. Higher polymer concentrations will stabilise the Taylor cone when the electric field pulls away the polymer jet. At lower polymer concentrations, a weaker Taylor cone forms and is subsequently interrupted resulting in beading when a voltage is applied. Besides the amount of polymer required to stabilise the Taylor cone, the flow rate used has to be sufficient to yield a continuous Taylor cone ( figure 5(B)).
The viscosity of the polymer solution is directly related to the MW of the polymer [28,33,46,47], structure of the polymer, additives (such as addition of salts) and solvent system used [20]. The effect on the solution viscosity must therefore be considered when changing the polymer or solvent system [22,48]. A previous study observed a correlation between polymer entanglement concentration (C e ; i.e. the extent of interaction between polymer chains) and bead formation [49]. Beading was generated at a minimum C , e 1 while uniform and defect-free fibres were observed at C . e 2 2.5 -Entanglement concentration increases with higher polymer concentrations. This explains the decrease in beading observed when increasing the polymer concentration.
The current limited knowledge of the effects of interacting electrospinning parameters on fibre morphology requires further understanding of the relationship between those parameters and bead formation. Understanding this relationship will make optimisation easier and less time consuming without having to perform many experiments where a single parameter is tested in turn. From the DoE model, it is evident that the concentration of the polymer solution significantly influences the fibre morphology whereas polymer concentration and flow rate both significantly influence the number of beads in the scaffold. A previous study on PLLA fibres [50] showed that only polymer concentration had a significant impact on the extent of bead formation which is in contrast to our findings reported here. This could be due to the lower levels of polymer concentrations explored by Patra and colleagues (4 and 7% w/v compared to 15 and 25% w/v in this study), a difference in the solvent system used (DCM:DMF compared to Ac:DMF in this study) and the size of the DoE of only 30 runs (compared to 81 runs in this study). We found that the mean fibre diameter was mostly influenced by the applied voltage and the interaction between polymer concentration and flow rate. These findings are not fully in accord with previous DoE studies. For example, Coles et al [51] observed a significant impact of polymer concentration on the resultant PLA average fibre diameter with potential and distance having much lower impact on the fibre diameters. However, their study excluded data from experiments which resulted in no fibre formation, thus reducing the range of morphologies considered by their model. Other studies using different polymers saw the same significant impact of polymer concentration on the fibre diameter with no significant impact of the voltage and flow rate/polymer concentration as reported here. Again, this could be due to selected data sets being considered in the prior studies rather than the whole sets of fibres, including the 'imperfect' ones. Differences in the previous work include the exclusion of fibre defect formation [52,53], using only visual observation as a measure of beading formation [54], not including flow rate in their DoE [54][55][56], not factoring in the difference in the characteristics of the polymers [52-54, 57, 58] and solvents used [56]. For example, Senthil et al looked at the impact of polymer concentration, applied voltage and flow rate (μl min −1 ) on the average fibre diameter of styrene-acrylonitrile copolymer (SAN) using DMF as solvent [58]. Clearly, even if using a similar solvent (DMF to Ac/DMF in this study), the electrospinning of a polymer such as SAN with very different physical properties to the PDLLA in our study, would greatly affect the viscosity, conductivity and surface tension of the electrospinning solution. Therefore, this model has limited application to other polymer/solvent systems, due to the impact of solvent differences (such as, volatily, permittivity) and polymer characteristics (such as molecular weight) within the electrospinning system [20,26,59,60]. However, this DoE model could be expanded with more parameters such as solvent systems and different isomers of PLA to investigate the impact further. 4D contour prediction plots were produced from the data generated in this study and from the observed interactions in the DoE model ( figure 6). These prediction plots show that uniform fibres or fibres with a very small number of defects could be obtained with mean fibre diameters between 600 and 1000 nm when using the highest polymer concentration (25% w/v; figure 6(A)). Fibre mats with a small number of beads per μm 2 were observed at the highest flow rate (5 ml h −1 ), at a PDLLA concentration of 25% w/v and a collector plate distance of 20 cm (figures 6(B) and (C)). Although this resulted in more defects, these parameters generated scaffolds with the lowest mean fibre diameters (650-700 nm; figure 6(i)). When defect-free fibres were apparent, the range of fibre diameters were found to increase to between 700 and 1000 nm ( figure 6(D)). An increased number of beads will occur at high flow rates and low polymer concentrations, with the highest number of beads observed at very low or high voltages and larger collector plate to needle distances (figures 6(B) and (C)). Higher mean fibre diameters will be observed when using low applied voltages. When defect-free fibres are required, a compromise must therefore be made as the mean fibre diameter will increase, whilst for smaller fibre diameters a small number of defects are likely to occur (figure 6(D), (i) and (ii)). It may therefore be appropriate to set tolerance limits for such beading defects with the knowledge that these will be unavoidable.

Conclusion
A DoE model is an effective approach to identify the most significantly impacting parameters and the effects of parameter interactions on resultant electrospinning PDLLA scaffold characteristics. Identifying these significant parameters can streamline optimisation experiments thus reducing the time spent on preliminary studies to find the required parameter combination for the desired fibres. By separating the morphological characteristics of the fibre scaffolds, it was clear that polymer concentration plays a vital role in influencing the final fibre morphology and extent of beading. Polymer concentration plays a more interacting role with the flow rate in terms of influencing fibre diameter. Applied voltage also had a significant effect on the fibre diameter, which confirms the importance of the formation of a stable polymer jet and Taylor cone during the electrospinning process. Matrix inversion on the regression equations (supplementary table S3, equations (S1)-(S6)) as well as interpretation of the 4D contour plots (figure 6) reported by the DoE can be used as guidance to identify the most appropriate parameter combination for the desired scaffold morphology and fibre diameter. For example, when uniform fibres are required, high polymer concentrations and low (1 ml h −1 ) to medium (3 ml h −1 ) flow rates are recommended. Production of fibre mats composed of small fibres of <700 nm in diameter will result in a small number of fibres exhibiting the 'beads-on-strings' morphology. In this way, these regression equations give the parameters needed to result in the desired fibre diameter and provide a valuable indication of expected bead formation, or avoidance, as well as morphology.