Electrospun silk-based nanofibrous scaffolds: fiber diameter and oxygen transfer

In this study, silk fibroin was extracted from cocoons of silkworms and fabricated into nonwoven mats by electrospinning method. A new model based on the group method of data handling (GMDH) and artificial neural network (ANN) was developed for estimation of the average diameter of electrospun silk fibroin nanofibers. In this regard, concentration, flow rate, voltage, distance, and speed of collector were used as input parameters and average diameter of the fibers was considered as output parameter. Two models were capable to estimate average diameter of fibers with good accuracy. The average absolute relative deviation for GMDH and ANN models was equal to 3.56 and 2.28 %, respectively. Furthermore, due to importance of oxygen delivery to site of injury to promote wound healing, continuity equation for mass transport was employed for prediction of oxygen profile in the system containing wound dressing and skin. The result showed that our prepared wound dressing is capable to pass the oxygen completely to the skin layer and is not acting as a barrier for oxygen delivery to wound site. Since average nanofibers diameter can influence the mat physical, mechanical and biological properties then this model may serve as a useful guide to obtain tailor made and uniform silk nanofibers at various combinations of process variables.


Introduction
Nanotechnology is a growing field of manufactured materials with sizes less than 1 lm, and it is particularly useful in the field of biology because these applications replicate components of a cell's in vivo environment. Nanofibers, which mimic collagen fibrils in the extracellular matrix (ECM), can be prepared from a host of natural and synthetic biomaterials and have multiple properties that may be beneficial to field of materials for medical applications. These properties include a large surface-area-to-volume ratio, high porosity, improved cell adherence, proliferation and migration, and controlled in vivo degradation rates. The large surface area of nanofiber mats allows for increased interaction with compounds and provides a mechanism for sustained release of antibiotics, analgesics, or growth factors into injured tissue; high porosity allows diffusion of nutrients and waste. Improved cell function on these scaffolds will promote healing. Controlled degradation rates of these scaffolds will promote scaffold absorption after its function is no longer required (Hromadka et al. 2008). All the mentioned properties can be influenced by fiber morphology and diameter (Min et al. 2004;Kim et al. 2003). It is mentioned in the literature that, those fiber morphology and diameter are depend on many parameters which can be divided into four main categories: polymer properties (molecular weight and solubility), polymer solution parameters (polymer concentration, solution viscosity, conductivity, surface tension, and etc.), processing conditions (applied voltage, nozzle-collector distance, feed rate, and needle diameter), and ambient parameters (temperature, atmosphere pressure, and relative humidity) (Nasouri et al. 2012).
Silk fibroin, derived from silkworm cocoons, has attracted the scientific community due to its good biocompatibility along with suitable mechanical property. Fibroin protein is the major constituent 75 % of the cocoon and the remaining 25 % is sericin protein Chirila et al. 2013). The molecular orientation makes this protein form a semicrystalline structure which contains two phases: highly crystalline antiparallel b-sheet structure and non-crystalline part (Wang et al. 2004). The crystalline part lead to increase the strength and toughness and the non-crystalline part contributes the flexibility and elasticity to the fiber . Over the past years, many studies have explored silk fibroin in various forms from its regenerated solution, including porous scaffolds, films and electrospun fibers (Wharram et al. 2010). The electrospun silk fibers with uniform micro or nano scale fibers have found utility in producing different biomaterials like wound dressing or scaffolds for a variety of biological applications such as bone, nerve and skin tissue (Valenzuela et al. 2012).
The lack of comprehensive and predictive models is sensing in the field of electrospinning. For instance, it can be characterized that when one of the parameters of electrospinning process is changed what will happen for diameter of nanofibers. The experimental observations show an increase or decrease in fiber diameter. But, it is not possible to predict the diameter magnitude without conducting experimental procedure. Consequently when a special diameter of fibers is required, it is necessary to change various parameters. However, experimental investigations are time consuming and almost expensive. In these situations, existence of predictive models will be so beneficial. Recently, studies have been carried out to determine the feasibility and to optimize the diameter of electrospun nanofibers with different type of mathematical relationships or models such as design of experiments or artificial neural network (Malallah and Nashawi 2005).
Artificial neural networks (ANN) are flexible modeling method which shown excellent performance for modeling different problems. Generally, ANNs are empirical mathematical tools which can model various data sets even in the cases that complex relation is existed between input and output parameters. This method provides flexible non-linear mathematical function mapping of a set of input variables into the output of network. Although ANNs are accurate systems for modeling different problems, but they have a disadvantage about their mathematical structure (Atashrouz et al. 2014a, b). Obtained mathematical structures for an ANN model are complicated and practical application of this type of models is not easy. To address this issue, group method of data handling type neural networks (GMDH-NN) was proposed. GMDHtype neural networks also known as polynomial neural networks has both accuracy in modeling and simplicity in mathematical structure (Atashrouz et al. 2014a, b;Dodangeh et al. 2014).
In this work, GMDH and ANN models are used for prediction of fiber diameter and also comparison of their performance as two different types of neural networks to prepare electrospun nanofibers from silk fibroin. In this regard, concentration, flow rate, voltage, distance, and speed of collector were used as input parameter and diameter of fibers was considered as output parameter. In addition, due to importance of oxygen delivery to site of injury and to promote wound healing process, it will be interesting to predict the oxygen profile within the electrospun mat and also the skin layer as an example. In this regard, a mathematical relationship based on a finite difference method was applied for predicting profile of oxygen in system containing electrospun mat and skin.

Preparation of silk fibroin
Silkworm cocoons were obtained from Gilan University (Rasht, Iran). First, cocoons of Bombyx mori were boiled for 45 min in an aqueous solution of 0.02 M sodium carbonate, and then rinsed with distilled water to completely remove sericin. After drying, the degummed silk fibers were dissolved in 9.3 M LiBr at 60°C for 4 h and then were dialyzed in a dialysis bag (12,000 MWCO, Sigma, USA) against distilled water for 3 days with changing water several times for the removal of the salt. Finally, the film of silk fibroin solution was prepared (Mirahmadi et al. 2013;Uttayarat et al. 2012).

Morphological observation by scanning electron microscopy (SEM)
The morphological investigation of the electrospun nanofibers was performed with scanning electron microscope (AIS2100, South Korea). The intact samples coated with gold for SEM witnessing. In the SEM photos, the fiber diameters were determined by means of Image J software and the results were given as the average diameter ± standard deviation.

FT-IR spectroscopy
FT-IR spectroscopy (Tensor 27, Brucker Optics, Germany) was used to determine the protein conformation in electrospun nanofibers of silk fibroin. The scan was collected under absorbance mode from 4000 to 400 cm -1 at 4 cm -1 scan resolution.

Group method of data handling (GMDH)
In the GMDH type neural network algorithm, by combination of more than two variables at a time a polynomial expression can be obtained for studied problem. The algorithm is tried to find most appropriate configuration of polynomials which best accuracy. The following series is the form of grand correlation multinomial in GMDH systems (Atashrouz et al. 2015a, b): where y 0 is the output, x i , x j , and x k are input parameters and l is the number of layers. In addition, a 0 , a ijk are coefficients of grand correlation multinomial which should be adjusted based on the algorithm. Now for an input vector of x = (x 1 , x 2 , …, x n ) the output of the multinomial is f i ðx 1 ; x 2 ; . . .; x n Þ. The parameters of Eq. 1 should be adjusted based on the real data. It is supposed that real data are show with r i . So, to optimize the grand correlation multinomial, square of difference between real data and calculated data of multinomial expression should be close together: where M is corresponds to the number of observations. Figure 1 shows a schematic of proposed GMDH model. For more detailed information about the group method of data handling, the reader is referred to the literatures (Atashrouz et al. 2015a, b).

Artificial numeral network (ANN)
Neural networks consist of arrays of simple active units linked by weighted connections. ANN consists of multiple layers of neurons arranged in such a way that each neuron in one layer is connected with each neuron in the next layer. The network used in this study is a multilayer feed forward neural network with a learning scheme of the back-propagation (BP) of errors and the Levenberg-Marquardt algorithm for the adjustment of the connecting weights. Neurons are the fundamental processing element of an ANN, which are arranged in layers that make up the global architecture. The ANN input is the first layer in the network through which the information is supplied. The number of neurons in the input layer depends on the network input parameters. Hidden layers connect the input and output layers. Hidden layers enrich the network for learning the relation between input and output data. In theory, ANN with only one hidden layer and enough neurons in the hidden layer, has the ability to learn any relation between the input and output data (Zhang et al. 1998). Transfer function is the mathematic function that determines the relation between neuron output and the network. The sigmoid transfer function is as follow: (3) is the number of inputs to the neuron. ''w i '' is the weight coefficient corresponding to the input ''x i '' and ''O pj '' is the output corresponding to the ''j'' neuron. For completion of this section, we illustrate the learning BP algorithm (Atashrouz et al. 2013).

ANN training algorithm
The Back-Propagation Algorithm is one of Least Mean Square methods, which is normally used in engineering. In a multilayer perceptron, each neuron of a layer is linked to all neurons of the previous layer. Figure 2 shows a perceptron with a hidden layer. Each layer output acts as the input to the next neurons. To train Multilayer Feed Forward Neural Network, Back-Propagation Law is used. In the first stage, all weights and biases are selected according to small random numbers. In the second stage, input vector X p = x 0 , x 1 , …, x n-1 and the target exit T p = t 0 , t 1 , …, t m-1 are given to the network, where the subscripts n and m are the numbers of input and output vector, respectively. In the third stage, the following quantitative values are calculated and transferred to the subsequent layer until it eventually reaches the exit layer.
The fourth stage begins from the exit layer, during which the weight coefficients are corrected.
where ''W ij (t)'' stands for the weight coefficients from node ''i'' to node ''j'' in time ''t'', ''g'' is the rate coefficient, ''d Pj '' refers to the corresponding error of input pattern ''P'' to the node ''j'' and ''O Pj '' is the output corresponding to the j neuron. ''d Pj '' is calculated by the following equations for exit layer and hidden layer, respectively (Browne 1997): Here, the R acts for k nodes on the subsequent layer after the node ''j''. In the learning process, there are several parameters that have influence on the ANN training. These parameters are the number of iterations, number of hidden layers and the number of hidden neurons. To find the best architecture of the model, best set of the aforementioned parameters based on minimizing the network output error should be chosen (Atashrouz et al. 2013(Atashrouz et al. , 2014a.
Prediction of oxygen profile in a system containing wound dressing and skin Figure 3 shows a schematic of considered problem for prediction of oxygen concentration. As can be observed, the system is consisted of two regions of wound dressing and skin. In the wound dressing region, oxygen is only diffused and equation of continuity takes following form: However, since in the skin layer due to existence of oxygen consumption by cells a reaction term should be added to equation: Input Layer

Hidden Layer
Output Layer The reaction rate kinetic function describing the overall oxygen consumption for metabolic is assumed to has the Michaelis-Menten kinetics. In this regard equation takes following form: where v m and k m are Michaelis-Menten parameters which v m is the maximum rate of oxygen consumption and k m is the concentration of oxygen when the rate of reaction is equal to 1 = 2 v m . The value of diffusion coefficient, specific oxygen consumption rate is equal to 2.54 9 10 -5 cm 2 s -1 and 3.33 9 10 -5 s -1 respectively (Lee et al. 2013;von Heimburg et al. 2005).
The oxygen concentration at the upper boundary of the wound dressing is equal to atmospheric oxygen concentration (C atm ) which has the following relation: The equal flux condition exists at the interface of the wound dressing and skin which is as below: where L 1 is the thickness of wound dressing. In addition, diffusion of oxygen to the right and left side of boundaries is equal to zero (rC ¼ 0).

Modeling of the fiber diameter
To develop of GMDH and ANN models, concentration, flow rate, voltage, distance, and speed of collector were used as input parameter and diameter of fibers was considered as output parameter. Variable parameters of electrospinning are tabulated in Table 1. Optimized GMDH model was obtained based on the experimental data and tabulated in Table 2. In addition, it should be noted that experimental data were divided into two sections: randomly, 70 % of experimental data were used for training the model and the rest 30 % were considered for testing the model. It should be emphasized that test data set is necessary because in the proposed model, maybe performance for estimation of train data would be satisfactory, but when the model is applied for a new sample the result is not accomplished with good accuracy. Figures 4 and 5 shows estimated data by GMDH and ANN models as experimental data. Both models indicate good performance and experimental data are keep close to diagonal line. The average absolute relative deviation for GMDH and ANN models was equal to 3.56 and 2.28 %, respectively. The average absolute relative deviation for model is defined as follows: Comparing the MARD% of both models indicates that ANN is more accurate than GMDH. In addition to the good accuracy for both models, it will be interesting to check the physical response of the models to different parameters. It means, as explained in previous section when all input parameters expect one of them are constant, the diameter of fibers can increase or decrease with changing the considered parameter. To check this possibility for GMDH and ANN models, the effect of various parameters on response of both models was considered. Figure 6 represents the effect of flow rate on the nanofiber diameter and the predicted results by GMDH and ANN models. As shown in this figure, both GMDH and ANN models have well prediction about increased diameter and desired behavior. However, the results of ANN model for higher flow rate of 0.8 cc/h are different from GMDH model. It should be noted that according to existing theories of neural networks, this type of model cannot be extrapolated and just used for the modeling and interpolation between the experimental data. This means that it is better to use the neural networks in the range of variables to reduce the possibility of incorrect predictions. Figure 7 shows variations of diameter as speed of collector. As can be seen GMDH model was predicted a descending trend with increasing the speed of collector, nanofiber diameter decreased, but ANN model shows different slop trend. Thus, it can be concluded, ANN model has a good ability to conformity on the experimental data, but physically this model is not able to predict the trend of decreasing the nanofiber diameter.
In addition Fig. 8 shows variations of diameter as collector distance. According to this figure, ANN model shows irrational prediction and with increasing distance, nanofiber diameter first decrease and then increase, but Output layer d ¼ À84:2446 þ 0:00342148X 5 À 24:1294X 2 þ 100:383X 4 2 þ 1:65187N 1 À 4:80264 Â 10 À6 N 3 1 X 1 , concentration; X 2 , flow rate; X 3 , voltage; X 4 , distance; X 5 speed of collector Fig. 4 Comparison between estimated data by GMDH model and experimental data Fig. 5 Comparison between estimated data by ANN model and experimental data GMDH model as well predict the behavior of experimental data, but have some deviations with them. This is due to the type of functions is intended to develop the GMDH model. In addition, another advantage of GMDH model in comparison to ANN model is its simple mathematical structure. While, ANNs have complex mathematical form which make their applicability difficult.
Experimental investigation of the electrospinning parameters on fiber diameter Figure 9 shows FTIR spectroscopy of regenerated silk fibroin and as expected from previous works peaks in the Amide I (1655 cm -1 ), Amide II (1538 cm -1 ) and Amide III (1239 cm -1 ) related to random coil structure of it are observed (Rousseau et al. 2004;Mirahmadi et al. 2013).
The SEM images of electrospun silk nanofibers under different concentrations are given in Fig. 10. It is seen that increasing in polymer concentration from 10 to 13 % results in fiber diameter increase ranging from 135 ± 16 nm to about 300 ± 29 nm. This can be explained by two reasons: first, increase in the amount of polymer in electrospinning jet and second, more interaction between polymer chains in solution which later is lead to more resistance of solution against pulling by electrical charges. This behavior is also observed in the study of Sukigara et al. (2004). Figure 11 shows SEM images of electrospun nanofibers of silk fibroin under different spinning distances. If other parameters were kept constant, increasing of spinning distance is lead to more evaporation of solvent and consequently decreases in diameter of fibers slightly. It should be noted that the morphology of electrospun nanofibers is not changed strongly. So spinning distance is not considered very effective parameter in the morphology of electrospun nanofibers (Sukigara et al. 2004;Zhou et al. 2009). Figure 12 shows the effect of increase in flow rate on fiber diameter. Obviously with increase in the flow rate, more volume of solution is exit from the needle. Therefore, diameter of the fibers is increased. The flow rate has great importance effect on the formation of Taylor's cone and with increase of flow rate is caused to uniform morphology of silk fibroin nanofibers. Consequently, flow rate is a key parameter for obtaining an appropriate structure in electrospinning (Megelski et al. 2002).

Oxygen profile in wound dressing and skin
Available analytical methods are not capable to solve the continuity equation when reaction rate of Michaelis-Menten kinetic is presented. Consequently, to solve equations and boundary conditions we should use numerical methods. In this regard, implicit finite difference method was considered for solving the problem and programming was conducted in MATLAB software. It should be noted that parameter of Michaelis-Menten kinetic was used from refs (Lee et al. 2013;von Heimburg et al. 2005).
The result of solving the equations is presented in Fig. 13. The figure is shown the oxygen profile after 1 day that wound dress is attached to the skin. As can be seen, oxygen can be suitably diffuse on the wound dress and reach to the skin layer for consumption by cells. It should be noted that last layer of skin is more sensitive to the oxygen because the magnitude of oxygen has its lower value in this region. The lower magnitude of oxygen in last layer of skin is logical because the problem is diffusion control. It means the only mechanism to transport oxygen is diffusion mass transfer which should be contest with the oxygen consumption of cells. However, the figure shows that diffusion of oxygen in our proposed wound dress is appropriate for transporting of oxygen to even the last layer of skin.

Conclusion
In the present study, electrospinning of silk fibroin-based nanofibrous mat under different conditions was performed and effect of various parameters on the diameter of nanofibers was investigated. For example, we observed that increasing the solution concentration from 10 to 13 % results in increasing fiber diameters from 135 ± 16 nm to about 300 ± 29 nm. Oppositely, increasing distance from 8 cm to 12 cm results in fiber diameters decreasing from 180 ± 31 to 158 ± 23 nm. Furthermore, two mathematical models based on the GMDH and ANN were developed for prediction of nanofibers diameter. Both models have good accuracy (MARD% = 3.56 % for GMDH and MARD% = 2.28 % for ANN model), ANN model has less error than GMDH model, but it has not good ability to predict the trend of increase and decrease of diameter with variation of electrospinning parameters. So GMDH model is more appropriate for modeling the similar research study. In addition, a mathematical model was proposed for prediction of oxygen profile in wound dressing which it can be concluded that the problem is diffusion control. In  addition, the rate of oxygen consumption by cells is more than transported equation with diffusion mechanism.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Fig. 13 Three-dimensional oxygen profile in wound dressing and skin system (relationship between oxygen concentration and skin ?wound dressing)