Correlating in vitro performance with physico-chemical characteristics of nanofibrous scaffolds for skin tissue engineering using supervised machine learning algorithms

The engineering of polymeric scaffolds for tissue regeneration has known a phenomenal growth during the past decades as materials scientists seek to understand cell biology and cell–material behaviour. Statistical methods are being applied to physico-chemical properties of polymeric scaffolds for tissue engineering (TE) to guide through the complexity of experimental conditions. We have attempted using experimental in vitro data and physico-chemical data of electrospun polymeric scaffolds, tested for skin TE, to model scaffold performance using machine learning (ML) approach. Fibre diameter, pore diameter, water contact angle and Young's modulus were used to find a correlation with 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay of L929 fibroblasts cells on the scaffolds after 7 days. Six supervised learning algorithms were trained on the data using Seaborn/Scikit-learn Python libraries. After hyperparameter tuning, random forest regression yielded the highest accuracy of 62.74%. The predictive model was also correlated with in vivo data. This is a first preliminary study on ML methods for the prediction of cell–material interactions on nanofibrous scaffolds.

The engineering of polymeric scaffolds for tissue regeneration has known a phenomenal growth during the past decades as materials scientists seek to understand cell biology and cellmaterial behaviour. Statistical methods are being applied to physico-chemical properties of polymeric scaffolds for tissue engineering (TE) to guide through the complexity of experimental conditions. We have attempted using experimental in vitro data and physico-chemical data of electrospun polymeric scaffolds, tested for skin TE, to model scaffold performance using machine learning (ML) approach. Fibre diameter, pore diameter, water contact angle and Young's modulus were used to find a correlation with 3-(4,5-dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide (MTT) assay of L929 fibroblasts cells on the scaffolds after 7 days. Six supervised learning algorithms were trained on the data using Seaborn/Scikit-learn Python libraries. After hyperparameter tuning, random forest regression yielded the highest accuracy of 62.74%. The predictive model was also correlated with in vivo data. This is a first preliminary study on ML methods for the prediction of cell-material interactions on nanofibrous scaffolds.

Introduction
Tissue engineering (TE) scaffolds are primarily engineered to: (i) stimulate cell-material interactions and adhesion, and extracellular matrix (ECM) deposition, (ii) allow adequate transport of biological factors to enable cell survival, proliferation and differentiation, and (iii) promote processes such as angiogenesis and reduce inflammation [1]. Scaffolds have been developed for the engineering of a range of tissues such as bone, skin, muscle, vascular, neural, cartilage and ligament [2]. Regardless of the tissue type, a number of key considerations are important when designing or determining the suitability of a scaffold for tissue regeneration, namely biocompatibility, biodegradability, sufficient mechanical strength to ensure integrity and interconnected porous structures and high porosity to allow cell penetration and tissue integration, vascularization for the transport of gases, diffusion of nutrient and growth factors [3][4][5].
Polymeric TE nanoscaffolds have experienced tremendous progress with a number of scaffolds being used in clinical setting. Materials scientists are now faced with more sophisticated challenges such as being able to match materials and scaffold design with wound healing requirements. For instance, a preliminary matching of in vitro data with scaffold's physico-chemical characteristics may offer a better comprehension of the cell behaviour and bring evidence-based data to scaffold design. The use of computational methods in three-dimensional printing techniques for scaffold design, fabrication and simulation has been studied [6]. To the best of our knowledge, no study in the literature, up to now, has specifically looked at modelling of cell-material behaviour on electrospun scaffolds. Scaffold material, structure and fabrication technique are the three main parameters which determine scaffold properties. Thus, designing a computational model which can integrate both material and structure data is a challenging task. It has been argued that machine learning (ML) can offer an indispensable tool and overcome challenges in the biomedical field involving complex heterogeneous data when conventional statistical tools have failed [7][8][9][10]. Predictive modelling is a probabilistic process that allows us to forecast outcomes, on the basis of predictors. The latter are features that come into play in the determination of the required result, i.e. the outcome of the model [11]. However, the huge potential of ML for predictive modelling in TE still remains mostly unexplored.
In this paper, we hypothesize that using in vitro cell culture data and ML approaches, we can reverse engineer scaffold performance and predict in vivo outcomes (scheme 1). We have used in vitro and physicochemical characterization data, namely number of cells (quantified during 3-(4,5-dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide (MTT) assay), fibre diameter, pore diameter, water contact angle and Young's modulus, generated from a series of experiments using biopolymer-based scaffolds, to build a predictive model for scaffold performance. Our goal is to predict the in vitro outcome, given the physico-chemical features. Six simplified regression models were used and their metrics were compared to determine the most accurate algorithm. The accuracy of the ML model was then assessed against in vivo experimental results.

Scaffold fabrication
Scaffolds were fabricated using the electrospinning method (bottom-up NE300 Laboratory scale electrospinner, Inovenso Company, Turkey). Parameters for electrospinning were varied depending on the polymers within the blend and on the blend composition so as to produce random bead-free fibres. PHB/KCG and PHBV/KCG fibres were produced as reported [12]. Electrospinning parameters for PSuc-based and cellulose-based fibres have been reported, respectively, by Chummun et al. [13] and Ramphul et al. [14]. The fabrication of PDX/KCG and PDX/FUC has been detailed in Goonoo et al. [15].
Most blend solutions were prepared by mixing two solutions (solution A and solution B) with the exception of PDX/PHBV, PDX/PSuc and PLLA/PSuc. Table 1 summarizes the mass of polymers as well as the volume of solvents used for the preparation of different blend solutions. PLLA/cellulose and PLLA/cellulose 1% nanosilica mats were obtained following the deacetylation of PLLA/CA and PLLA/CA 1% nanosilica mats. Briefly, the acetylated mats were immersed in 0.05 M NaOH solution for 48 h at room temperature followed by washing with distilled water. The average fibre diameters of the electrospun mats were determined by scanning electron microscope (SEM) as reported previously [13,14,16]. Average pore diameters were determined by measuring the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201293 diameter of at least 50 different pores using ImageJ software. Pores were identified as areas of void space bounded by fibres on all sides at or near the same depth of field. The shortest diagonal axes were measured and averaged together to obtain pore diameter. SEM images were taken using a Tescan Vega 3 LMU microscope (CBBR, Mauritius) with accelerating voltage of 30 kV.

Mechanical properties
A Universal Instron Tester 3344 (Instron, USA) was used to measure the tensile properties of the electrospun mats at 25°C. The electrospun mats cut in rectangular shapes (4 × 1 cm) were clamped with gauge length set at 1 cm and strained at a rate of 10 mm min −1 using a load cell of 2 kN.

Water contact angle
The contact angles of the electrospun mats were measured using Milli-Q water as a probe liquid with a Krüss Drop Shape Analyzer DSA 25 (Advanced Lab GmbH, Germany). Static contact angle data based on the sessile drop method were acquired immediately after deposition of a 2 µl drop on at least three positions for each sample and are stated as the arithmetic mean.

In vitro evaluation by MTT assay
The MTT assay was conducted using L929 mouse fibroblasts (Sigma-Aldrich, ECACC certified) after 7 days as reported earlier [12] and absorbance values were measured at a wavelength of 540 nm using a Thermo Scientific Varioskan LUX Multimode Microplate Reader.

In vivo biocompatibility tests
Surgical procedures were approved by the Animal Ethics Committee of CYROI, La Réunion, France (APAFIS Number: 2018052311219125V3). Wistar albino rats (8-10 weeks old), both males (308-407 g) and females (189-252 g), were used for in vivo studies (n = 5). Prior to implantation, the scaffolds were disinfected in 70% ethanol solution overnight. The rats were anaesthetized by isoflurane inhalation and the dorsal region was shaved followed by 1-2 cm dorsal incisions (two on each side of the spine). All the scaffolds were implanted as stacks of three each with dimensions 0.5 × 1 cm. The rats were provided with regular supply of food and water throughout the study period. They were also monitored daily for any signs of inflammation. The scaffolds were removed and the neighbouring skin tissues were harvested for histology analysis after two and four weeks, respectively. Histological analysis of tissues was performed using Masson's Trichrome staining. The stained slides were examined using light microscopy (Nanozoomer 560 digital slide scanner, Hamamatsu). For quantitative analysis of angiogenesis, the number of blood vessels were counted within a fixed distance of 100 µm from the explanted scaffold and the results were expressed as the number of blood vessels per mm 2 for each scaffold [17].

Data preparation and supervised machine learning procedure
Before applying statistical learning to scaffold data, standard pre-processing was performed. Data were collected, cleaned and rearranged in a format appropriate for meaningful statistical analyses. The dataset used for this study contained 182 observations, four numeric features characterizing the scaffolds-fibre diameter, pore diameter, water contact angle and Young's modulus-and the target variable: number of cells. The variables were scaled by the MinMaxScaler method in order to bring the data on a relatively similar scale and close to the normal distribution. Low variance filter and high correlation filter using correlation matrix with Pearson correlation were applied on the features. Feature selection was performed with random forest regressor, recursive feature elimination (RFE) and forward feature selection (FFS) [18].
For predictive modelling settings, a data matrix denoted X is needed, and y as the target variable (to predict). We used 80% of the data for training, and the remaining 20% was used for testing. Six regression algorithms, namely linear regression, support vector regression (SVR), random forest regression, lasso regression, decision tree regression and k-nearest neighbour (k-NN) regression were applied to the data and compared to obtain the best performance. The six ML algorithms have undergone a training process to find patterns in the training data that map the input parameters to the target. Each algorithm learned from the training dataset containing both inputs and outputs to generate a ML royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201293 model. Each ML model was then tested using the test dataset. This metric has allowed us to evaluate how the model might perform against unseen data. To optimize the training phase, hyperparameter tuning was performed for each model in order to improve its accuracy. All codes were implemented in Python 3 using Seaborn/Scikit-learn libraries [19].

Importing Python libraries and loading the dataset into a data frame
Required libraries numpy and pandas were imported, and the read_csv method of pandas was used to load the dataset 'scaffold_data.csv' into a pandas data frame df, as shown below.

Splitting the dataset into training and test sets
Before feeding data to a regression model, some pre-processing is needed. A variable X was created from the initial dataset df to contain the four scaffold parameters (i.e. features): fibre diameter, pore diameter, water contact angle and Young's modulus. A variable y was created to represent the target variable, number of cells. The train_test_split function of Scikit-learn was used to split the data into training and test sets. test_size = 0.2 indicates that 20% of the initial dataset was used for testing, and 80% for training. random_state ensures reproducibility and X_train, X_test, y_train and y_test were defined for the ouput of train_test_split. from sklearn.model_selection import train_test_split X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, As variance is range dependent, feature scaling is required on the variable X before training the data. The MinMaxScaler method was used to scale and translate each feature individually such that it is in the given range on the training set, i.e. between zero and one. from sklearn.preprocessing import MinMaxScaler scaler = MinMaxScaler() scaler.fit(X) X_train = scaler.transform(X_train) X_test = scaler.transform(X_test) 2.6.3. Creating a random forest regression model and fitting it to the training data X_train and y_train obtained above were used to train a random forest model from RandomForestRegressor. The fit method was used to pass the parameters defined below. The output of this step describes a large number of parameters for the random forest model. These parameters are tuneable, and can be adjusted to optimize the accuracy of the model.

. Prediction
Once the model is trained, it is ready to make predictions. The predict method was used on the model and X_test was passed as a parameter to get the output as y_pred.

Statistical analysis
All group data were reported as the arithmetic mean ± standard error of mean. Statistical analysis was performed with GraphPad Prism software. Statistical differences were analysed using unpaired Student t-test. Differences were considered statistically significant for p < 0.05 ( Ã significant at p < 0.05, ÃÃ significant at p < 0.01, ÃÃÃ significant at p < 0.001, ÃÃÃÃ significant at p < 0.0001). Regression metrics such as accuracy measure, Spearman correlation and mean absolute percentage error (MAPE) were computed to evaluate the predictive accuracy of each ML model. Correlation between the actual and predicted output data was performed to assess the adequacy of each ML model.

Results and discussion
3.1. Physico-chemical characterization of nanoscaffolds: fibre diameter, pore diameter, water contact angle and Young's modulus In this study, we have used nanofibrous scaffolds generated through the electrospinning process. We have previously reported on the synthesis of these scaffolds which have been assessed for skin TE both in vitro and in vivo [12][13][14][15]17]. These scaffolds are engineered using blend solutions of various biopolymer combinations (table 2) and we have shown that the cell-scaffold response depended on the combination of materials as well as the physico-chemical properties of the scaffolds namely pore diameter, fibre diameter, water contact angle and Young's modulus (table 2). Fibre diameter and pore diameter were determined using ImageJ software. The data are presented as arithmetic mean ± standard error of mean. The different combinations of biomaterials used in this study will allow the testing of the robustness of the ML models to relate biological data with the physico-chemical characteristics of the scaffolds. L929 fibroblasts cells were grown on the scaffolds and cell proliferation was assessed using MTT after 7 days. Fibroblasts were chosen due to their crucial role in tissue regeneration in general across all tissue types. They are responsible for the proliferative phase in skin TE and are the major producer of ECM for the growth of keratinocytes and endothelial cells. Fibroblasts and macrophages are the first cells which interact with the surface of TE scaffolds in wounds and thus, determine the fate of the scaffolds as an efficient medical device.
Our experimental work indicated that there is no linear relationship between physico-chemical parameters and in vitro results. For instance, an increase in pore diameter did not translate directly to an increase in cell number but would influence contact angle and mechanical property, and together these will regulate cell-material interaction. The right balance between these properties needs to be determined to predict scaffold performance in vitro. Another observation that we made concerned the fibroblast morphology and behaviour (figure 1). For instance, changing the biopolymer from KCG to PSuc to cellulose caused a change in L929 morphology from elongated to dendritic to flat dense clusters, respectively. Each morphology change is an indication of a different type of cell-material interaction. Thus, feature selection was performed to select the predominant physico-chemical parameter which will impact not only cell behaviour but also the other physico-chemical parameters.
The main feature of ML is learning from experience. In the first step of supervised learning, a labelled training input dataset is fed to a ML algorithm. With the training dataset, the system adjusts itself by modifying parameters to create a logical model. The performance of the built model is then tested with a test dataset which the model has never seen before and the accuracy of the ML model is evaluated. Six supervised ML algorithms were used to predict scaffold performance with in vitro cell culture data.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201293 [ R d and the output variable y takes continuous values ðy [ RÞ. The goal of supervised learning algorithms is to analyse the training data and produce a function in order to predict the correct outcome for newly presented input data. The function used to connect input values/features x to a predicted output y, is created by the ML model during the training phase. The obtained function is then tested by how accurately it performs on unseen data assumed to follow the same distribution as the training data [20].
Data exploration and feature selection were performed on the 13 scaffold families (table 2), representing a dataset of 182 observations and four features ( pore diameter, fibre diameter, water contact angle and Young's modulus). Data were collected and regrouped to create an initial dataset df where each row is an observation (scaffold) and each column represents one feature (scaffold parameter) (electronic supplementary material, table S1). The initial dataset was loaded as a data frame and a variable X was created, including the four scaffold parameters (fibre diameter, pore diameter, water contact angle and Young's modulus). A variable y was also created and represented the target variable number of cells.

Low variance and high correlation filters
In ML, variances are calculated to make generalizations about a dataset, aiding in the understanding of data distribution. The variance of a variable is a measure of the average variation of values in the distribution with respect to the mean, given by where σ 2 is the variance, N is the total number of observations, X is the value of an individual observation and µ is the mean. Low variance filter calculates the variance of each variable in our dataset and removes low variances (below a given threshold, set at 10 −3 in our case) that would be of no significant contribution to the ML model. Variance calculations with scaled data was 0.017 for fibre diameter, 0.017 for pore diameter, 0.023 for Young's modulus and 0.149 for water contact angle.
High correlation filter calculates the correlation between independent scaled numerical variables. If the correlation between two variables is high (value closer to 1), it indicates that the variables have similar trends and are likely to carry similar information which can bring down the performance of some ML models. Figure 2 shows the Pearson correlation matrix performed on the four features (fibre diameter, pore diameter, water contact angle and Young's modulus). Coefficient correlations were between −0.19 and 0.42, i.e. no variable was highly correlated to another one.

Feature selection with random forest
Random forest is one of the most widely used algorithms for feature selection and allows us to compute the importance of each variable on the decision tree. The latter is a set of tests that are hierarchically organized and consists of nodes for testing features, edges to represent the outcome of each test and leaf/terminal nodes indicating the decision taken after computing all features. Random forest algorithm generates a large set of randomized decision trees to predict the target, and each feature's usage statistics are used to find the most informative subset of features in the dataset. Every decision tree has high variance, but when combined together and run in parallel, the resultant variance is low as random forests train each individual decision tree independently on different bootstrapped samples of the training data and hence, the output does not depend on one decision tree but multiple decision trees.
In the case of a regression problem, the predictions of each decision tree are averaged to make an overall prediction as final output (figure 3a). This technique is commonly known as bagging.
A score was calculated for each feature and the most predictive features were those with the highest scores. The feature importance graph showed that fibre diameter was the most important feature (0.39) while water contact angle and Young's modulus were the least important features (0.14 and 0.1, respectively) (figure 3b). This indicated that fibre diameter played an important role in dictating cell growth. The fibre surface is the first point of contact of cells as they first attach themselves to the fibre and exert the push-pull effect to be able to move and proliferate. Pore diameter, the second important feature, is responsible for cell penetration in scaffolds as well as the circulation of nutrients for cell growth. It is important to note that cells experience mechanical properties at the nanoscale and microscale, whereas mechanical properties of scaffolds are measured on the macroscale. Fibroblast cells are known to attach themselves on surfaces with water contact angle in the range of 60-80° [21]. This may explain the least impact exerted by this parameter as the scaffolds exhibited water contact angles in the range of 25-141.3°. The influence of this feature may vary with another cell type.

Recursive feature elimination and forward feature selection
RFE and FFS are applied to datasets with small number of input variables as is the case in this present study. RFE uses an accuracy metric to rank the features according to their importance. Linear regression was selected as model with four features and RFE gave '1' and 'True' for all four features. Three optimum features were found by the RFE method: fibre diameter, pore diameter and Young's modulus (score = 0.063). FFS is the opposite process of RFE as it tries to find the best features which can improve the performance of the model. In this study, the variables fibre diameter ( p = 0.003) and pore diameter ( p = 0.004) were retained as the most important parameters. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201293 RFE and FFS scores confirmed the scores computed with random forest regressor. Among the three feature selection methods performed on the data, fibre diameter and pore diameter were identified as the most relevant features in the dataset. These two parameters could contribute most to the construction of our predictive model.

Evaluation and selection of the predictive model
After the training process, model testing is performed to determine the accuracy of each ML model. In general, for a selected ML model, the scores generated on the different test sets are not enough to establish a quantitative score to select and assess the robustness of the model. Hyperparameters are predefined parameters fixed before the training process. They can be tuned and optimized to achieve the maximal performance of a model. Thus, model selection is performed by maximizing the crossvalidation score. Hyperparameters are tweaked in a principled way and a selection between different models is available.
Six regression models, namely linear regression, SVR, lasso regression, random forest regression, decision tree regression and k-NN regression, were employed to predict scaffold performance with in vitro cell culture data. Each ML algorithm was trained on an 80% random split of the dataset, while royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201293 the remaining 20% of the data was used to test the cross-validated model and evaluate the accuracy of the ML algorithm. For the ML models to provide optimal prediction performance, hyperparameter tuning was performed with the GridSearchCV estimator. For each ML model, cross-validation scores were computed for all hyperparameter combinations to find the best one. While selecting a regression model, computing metrics like accuracy measures and error rates are important to evaluate the prediction accuracy of the model. The MAPE is the most common statistical measure used to forecast error and is calculated as the average of the absolute percentage errors of predictions. MAPE can be formalized by the following mathematical expression where n is the size of the sample, P t is the value predicted by the model for time point t and A t is the value observed at time point t. The measure of the accuracy of the model is then calculated as accuracy = 100 − MAPE. Table 3 shows the accuracy score, MAPE and Spearman correlation coefficient to the corresponding regression method after hyperparameter tuning. Random forest regression gave the best performance with a decent accuracy of 62.74% and a high degree of Spearman correlation (0.64) implying that the actual and predicted values had similar directional movement, i.e. when the actual values increased, the predicted values also increased. Compared to the other models, the random forest method has shown the lowest MAPE (37.26%), indicating a reasonable prediction of our model [22]. The following set of hyperparameters were adjusted for the random forest method: n_estimators (number of trees in the forests, set at 200), max_features (number of features to consider for splitting a node, set to sqrt(n_features)), max_depth (maximum number of levels in each decision tree, set at 40), min_samples_split (minimum number of samples required to split an internal node, set at 2), min_samples_leaf (minimum number of samples required to be at a leaf node, set at 4) and bootstrap (using bootstrap samples when building trees, set to True). The five other ML models yielded accuracy scores lower than 60% indicating that these statistical models were not presenting accurate mapping between the inputs and predicted output (table 3).
The correlation between the actual and predicted output data is another important factor when designing a predictive model. We evaluated the adequacy of the model by plotting the actual values from the test dataset against predicted values from the random forest model ( figure 4). The plot displayed a statistically significant decent fit ( p < 3.7 × 10 −4 ), inferring a potential relationship between the predictors and the outcome that can be improved with more values closer to the fitted regression line. R-squared (R 2 ) represents the percentage of variance explained by covariates in the model and measures the proportion of variability in the outcome that has been explained by the model. It is ranged between zero and one. Usually, the higher the value, the better the model is able to explain the variability in the outcome, depending on the quality of the data. In our case, R 2 (0.3073) showed that more than 30% of the variance in the data is explained by the model. It is hypothesized that a better prediction accuracy will be seen by increasing the number of observations in the training set, and having more features in our dataset.
Fibre diameter values between 0.13 and 1.35 µm (including standard error of mean) were on/closest to the fitted line. An analysis of the experimental data indicated that these fibre diameters corresponded to the presence of biopolymers in the ratio range 10-30% in the polymeric blends with PDX, PHBV, PHB or PLLA (electronic supplementary material, table S2). 3.4. Preliminary correlation with in vivo data: impact of fibre diameter on angiogenesis The ML findings were correlated with in vivo biocompatibility experiments conducted on Wistar rats. Two families of scaffolds (PDX/CA and PLLA/cellulose) were implanted in the dorsal region of the rats and after two/four weeks, the scaffolds and surrounding tissues were removed and analysed. Angiogenesis is an important phenomenon in the wound healing process. The number of blood vessels was counted within a fixed distance of 100 µm from the explanted scaffold in a measured surface area (electronic supplementary material, figure S1). In this quantification study, a blood vessel was identified as one with a well-defined lumen consisting of red blood cells and an intact wall. The results showed a direct correlation between fibre diameter and angiogenesis whereby a decrease in fibre diameter resulted in an increase in number of blood vessels adjacent to the scaffolds, indicating enhanced angiogenesis (table 4).

Conclusion
We have shown for the first time the implications of incorporating ML methods to nanoscaffolds performance. A simplified predictive model was constructed and six regression algorithms were tested using the Seaborn/Scikit-learn Python toolkit. Hyperparameter tuning was performed to choose a set of optimal hyperparameters for each model in view of obtaining the best performance of each ML model, and random forest regression was found to be the model with the highest accuracy. Fibre diameter and pore diameter emerged as the two physico-chemical parameters which impacted more on the MTT values. Algorithms from transfer learning and reinforcement learning for increased accuracy are now being considered to increase the robustness of our model in predicting in vitro outcome. The final outcome would be a generalizable model that uses physico-chemical parameters of scaffolds as input and generates predictions of cells proliferation during in vitro cell culture in new scaffolds. The impact of different cell types on the model will also be studied.