An artificial-vision- and statistical-learning-based method for studying the biodegradation of type I collagen scaffolds in bone regeneration systems

Texture analysis

The approach proposed in this study is to quantitatively determine cell growth and differentiation and the degree of degradation of type I collagen due to a biological activity. The first step for achieving these objectives is to associate a numerical vector of features with each pixel. In fact, once the different ROIs were identified by expert personnel, readers can see that the different regions can be differentiated according to the texture of the image. For this purpose, scale-invariant texture filters were applied to 2D images by means of the ImageJ software. Image texture can be defined as a spatial arrangement of the color or intensities of an image in a given region, whereas texture analysis is often performed to separate different regions of an image. Therefore, in this paper, we propose to apply an image texture analysis consisting of extraction of representative features, such as a structure tensor and entropy, to identify different ROIs.

Once the characteristic vectors are extracted, they will serve as a training sample for the classifier, in this case, the random forest algorithm.

1. Structure tensor: This is a matrix representation of the image’s partial derivatives defined as second-order symmetric positive matrix J: (1)${\displaystyle J=\left[\begin{array}{cc}\hfill 〈f_x,f_x〉w\hfill & \hfill 〈f_x,f_y〉w\hfill \\ \hfill 〈f_x,f_y〉w\hfill & \hfill 〈f_y,f_y〉w\hfill \end{array}\right]}$ where f_x and f_y are images of the partial spatial derivatives: $\frac{\partial f}{\partial x}$ and $\frac{\partial f}{\partial y}$, respectively (Budde & Frank, 2012). From this matrix, all major and minor eigenvalues are separated for each pixel and channel in the image: (2)${\displaystyle T\left(v\right)=\lambda v}$ where λ is the eigen values, and v denotes the eigen vectors.

2. Entropy: For its calculation, a circle of radius r is drawn around each pixel. The image intensity histogram corresponding to the separated circle is obtained as binarized image fragments. Finally, entropy is calculated for each particle, where p is the probability of each chunk in the histogram corresponding to each channel of the image, both in RGB and in Hue, Saturation Brightness formats.

The random forest classifier for identifying different ROIs

It is possible to apply a wide variety of supervised classification models based on the descriptors that were extracted from the images in order to identify different ROIs. For the specific case of this research, a random-forest–type classifier (Criminisi, Shotton & Konukoglu, 2012) was chosen. It has the advantages of operating on attributes (or modular characteristics), of preventing overfitting of certain classes, and an optimized computational cost.

In recent years, decision trees proved to be some of the most promising techniques in machine learning, computer vision, and medical analysis (Criminisi, Shotton & Konukoglu, 2012). Random forest classifiers operate by constructing several decision trees (predictive processes that map observations of an item to conclusions about the objective value of the item) in the training phase, to then produce a class fashion (by its nature as a classifier) for each tree.

Image segmentation

Now, the benefits of the technology based on artificial vision allow us to determine the biological behavior of cells such as cell growth and differentiation as well as degradation of biomaterials such as type I collagen. Image segmentation is a reference within the support techniques for achieving this goal. In general, image segmentation is defined as the process of dividing an image into different segments or groups of pixels that share certain common characteristics (Mallik et al., 2011; Janeiro-Arocas et al., 2016). This is a fundamental task in image analysis for the detection of objects (Mallik et al., 2011; Janeiro-Arocas et al., 2016). In fact, it is a central problem in many studies dealing with image analysis (Meijering, 2012). It is important to note that the type and quality of the acquired images influence the success of cell segmentation (identification and separation of objects) (Kasprowicz, Suman & OToole, 2017). Therefore , adequate quality of images is one of the existing requirements of this set of techniques. In fact, an important feature of image segmentation is that although it is conceptually simple, it lacks generality; consequently, it cannot be implemented reliably and easily for all cell lines, image modalities, and cell densities without preprocessing of images (Alanazi et al., 2017).

Limiting factors for the image segmentation process are, e.g., the type of objects to be analyzed in the image, the aims pursued by the research team, and limitations of the knowledge of the technician in charge of the segmentation process. All these factors necessitate devising specific strategies for the processing and analysis of images for each particular problem. Despite the absence of a universal segmentation procedure (Alanazi et al., 2017), today it is possible to analyze 2D images of the behavior of stem cells in vivo, with information on sequential growth and differentiation as a function of time, via a set of different segmentation techniques. Among the most popular segmentation and image-processing techniques are thresholding, region-growing methods, edge detection, and Markov random field color algorithms (Grys et al., 2017). In this work, we propose a segmentation method based on a random forest classifier trained on vectors of attributes and on properties obtained from filters that are applied to labeled ROIs in the images. In fact, different studies have revealed the great potential of the random forest method as a support tool in the segmentation process. Thus, Schroff, Criminisi & Zisserman (2008) demonstrated the potential of random forest as a technique for segmentation processes of objects present in images. They worked with the MSRC image dataset, incorporating both features that describe objects locally and those that characterize the environment of those objects. In the medical domain, Dhungel, Carneiro & Bradley (2015) identified formation of unusual masses in mammograms. These masses varied in size, texture, and the area that delimits them, merging with the rest of the tissues of the breast. The classifiers based on deep belief networks (first filtering in which suspicious areas are located) and random forest (last filtering) were applied in a cascade. This approach significantly reduced the number of false positives. Furthermore, Khan, Hanbury & Stoettinger (2010) proved that the random forest technique is highly effective at segmenting human skin images. Those authors used random forest in combination with Improved Hue, Luminance and Saturation color space and compared this classifier with others based on Bayesian networks, multilayer perceptron, support vector machine, AdaBoost, naïve Bayes, and Radial Basis Function neural networks. The classifier surpassed all other techniques, yielding 0.877, 0.738, and 0.740 in terms of accuracy, precision, and recall, respectively. Finally, Ghose et al. (2012) demonstrated the effectiveness of segmenting magnetic resonance images of human prostates. To this end, those authors propose to use a framework of decision trees in order to obtain a probabilistic representation of those voxels that define the prostate.

The supervised classification process was applied by means of the Weka software (Hall et al., 2009). This computational tool allows for the application of a wide range of techniques from the field of data mining, e.g., data visualization techniques, feature selection, unsupervised classification or clustering algorithms, supervised classification such as random forest, and regression models (Frank et al., 2004).

Regression model fitting

Once the vectors of features related to the growth and cellular differentiation of osteocytes were obtained through the process of image segmentation (mentioned above), the degree of degradation of type I collagen was estimated via application of statistical regression models. In addition, it is important to note that the extracellular matrix ROI is extracted from the image and its features. They are the area, objects’ circularity (4πarea∕perimeter²), and the intensity histogram mode, among others that are modeled as a function of time to estimate cellular growth. The relation between collagen mass loss and the area or circularity was studied here to model the degree of degradation of collagen as well as the biological activity.

Two main regression approaches were employed: nonparametric and parametric nonlinear regression. In the text below, the two aforementioned approaches are briefly introduced.

Nonparametric methods have attracted the attention of academia and industry (Hayfield & Racine, 2008). They have relevant advantages, e.g., they do not assume any parametric function that constrains the relation between variables or even a specific distribution of those variables of interest (Tarrío-Saavedra et al., 2011). Indeed, nonparametric methods provide estimates of flexible models that do not impose any prespecified function (Maity, 2017), thereby allowing investigators to model a wide range of possible complex nonlinear functions (Shokrzadeh, Jozani & Bibeau, 2014). There are many types of nonparametric models dedicated to regression tasks, for example, kernel smoothing (Wand & Jones, 1994), local polynomial regression (Breidt & Opsomer, 2000), regression splines (Wood, 2006; Mammen & Van de Geer, 1997), and LOWESS (Cleveland & Devlin, 1988).

Specifically, penalized regression spline modeling within the framework of nonparametric GAMs (Wood, 2006) is proposed here to obtain information about the type of dependence between collagen mass loss and time. For the sake of simplicity, a penalized B-spline basis was fitted. A generic GAM as a function of two linear predictor variables (X₁ and X₂) and two smooth predictors (T₁ and T₂) can be defined as follows: (3)${\displaystyle Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+s_1\left(T_1\right)+s_2\left(T_2\right)+s_{12}\left(T_1,T_2\right)+ϵ}$

where the effects on the response Y of X₁ and X₂ are linear, whereas the effects of T₁ and T₂ are only assumed to be smooth. This model additionally allowed us to include linear or smooth effects for the interaction between variables, as is the case for S(T₁, T₂). The discrepancy between Y and Y estimates, $\hat{Y}$, was measured by the model residuals, ϵ. When the Y variable is discrete, e.g., binomially or Poisson distributed, Y is replaced by g(η). In the present case, Y = MassLoss and X = Time; thus, (4)${\displaystyle \hat{Y}=E\left(Massloss|Time\right)={\beta }_0+s_{Time}\left(Time\right)}$

As mentioned above, package mgcv provides computational tools to fit GAMs. This library uses penalized regression splines for the estimation of smooth effect s of each predictor variable x. To define the concept of splines, consider the set of nodes a < t₁ < t₂ < ⋯ < t_L < b in interval [a, b]. A spline function is a polynomial function in each subinterval of a certain degree (e.g., three in the case of a cubic spline), which is continuous in the knots up to order degree L − 1. With the aim of determining a spline function, we must know L + degree + 1 values, known as the degrees of freedom of the spline family defined for the set of knots. The degree + 1 number is known as spline order.

Thus, we can write the smooth effect as (5)${\displaystyle s_j\left(x\right)=\sum _{k=1}^{L+degree}{\beta }_{jk}{\varphi }_{jk}\left(x\right)={\beta }_j^{\prime }{\varphi }_j}$ for the basis composed of L + degree functions ϕ_jk that define the B-splines, where L is the number of knots (10 by default) and degree = 3 (cubic regression splines). Given the s_j definition, $\hat{Y}=E\left(Massloss|Time\right)$ is now a parametric mean function with parameters $\beta ={\left({\beta }_0,{\beta }_{Time}^{\prime }\right)}^{\prime }$ and predictor Time that define the intercept and the splines that define s_j.

The penalized least squares objective function for estimating β is (6)${\displaystyle \parallel Y-X\beta {\parallel }^2+\sum _j{\lambda }_j{\beta }_j^{\prime }{B}_j{\beta }_j}$ The values of λ_j are selected by an iterative algorithm to minimize a generalized cross-validation criterion.

Parametric models provide the transfer function that determines the dependence relation between critical variables of a process. In addition, the constitutive parameters usually have a physical meaning. Thus, the chemical, physical, and biological processes of materials and substances can be characterized and compared by studying the parameter values.

To illustrate a parametric-model expression in a simple way, we can define (7)${\displaystyle Y=m\left(t_i\right)+{\epsilon }_i,\mathrm{with}i=1,2,\dots ,n,\mathrm{and}\mathrm{E}\left({\epsilon }_i\right)=0}$

where Y is the response variable to be modeled, and m is the regression function fitted to Y. The latter can be the mass loss, area, and circularity, among other features critical for collagen degradation and cell growth. Finally, ε_i are the independent, identically distributed model residuals. In a fixed design, 0 ≤t₁ ≤t₂ ≤… ≤t_n, the parametric functions can be defined as $\mathcal{M}=\left\{{\mathrm{m}}_{\mathrm{\theta }}∕\mathrm{\theta }\in ⊝\right\}$, where θ is the parameter vector that defines the regression model, and ⊝ is a subset of ℝ^k. The linear models are the simplest and most typical parametric models, but in chemistry and biology domains, nonlinear models are also important and useful. In this regard, there are many processes characterized by exponential and sigmoid types of relations between critical variables such as material degradation (Robles-Bykbaev et al., 2018; Ríos-Fachal et al., 2014), crystallization (López-Beceiro, Gracia-Fernández & Artiaga, 2013), oxidation (Tarrío-Saavedra et al., 2013), and material dimensional variation due to magnetic changes (Tarrío-Saavedra et al., 2017). In fact, the requirements of biological growth have popularized sigmoid monotonic nonlinear functions, characterized by exponential growth and a point of inflection from which the speed of growth gradually decreases until reaching a saturation zone (Kahm et al., 2010; Román-Román & Torres-Ruiz, 2012; Kingsland, 1985; Tarrío-Saavedra et al., 2014). Some of the most popular functions are those based on logistic and Gompertz functions:

1. Gompertz function: (8)${\displaystyle y\left(t\right)=Aexp\left[-exp\left(\frac{\mu e}{A}\left(\mathrm{\lambda }-t\right)+1\right)\right]}$

2. Logistic function: (9)${\displaystyle y\left(t\right)=\frac{A}{1+exp\left(\frac{\mu -t}{scal}\right)}}$

where A is the asymptote, µ represents the inflection point or time corresponding to the maximum rate of change, λ is the delay until the beginning of exponential growth, and scal or the scale parameter accounts for the growth (or decay) rate of the function. The logistic function is used to explain a great number of real processes related to biology, chemistry, physics, informetrics, and economy, among other fields, taking into account that it is their underlying model (Limoges, 1987; Román-Román & Torres-Ruiz, 2012). Furthermore, nonlinear mixture models are created by the addition of two or more nonlinear functions. They serve for identifying, separating, and estimating overlapping processes often accompanying further kinetic analysis (Francisco-Fernández et al., 2012; Sánchez-Jiménez et al., 2013; López-Beceiro et al., 2010; López-Beceiro et al., 2012). This is the case for the mixture logistic model: (10)${\displaystyle m\left(\theta ;t\right)=\hat{Y}\left(t\right)=\sum _{i=1}^k\frac{A}{1+exp\left(\frac{{\mathrm{\mu }}_i-t}{{scal}_i}\right)},}$

where k denotes the number of logistic components, t is usually time or even temperature, and θ is the vector of parameters θ = (A_i, μ_i, scal_i). Logistic mixture models have been successfully applied in many studies on overlapping processes (e.g., degradation and crystallization processes) of different materials such as wood, polymers, composites, ceramics, and metal-organic frameworks (Rios-Fachal et al., 2013; López-Beceiro et al., 2011; López-Beceiro et al., 2010; Pato-Doldán et al., 2012; Francisco-Fernández et al., 2012).

Exploratory correlation analysis

Identifying and characterizing dependence relations between features extracted from an image, time, and the degree of degradation (collagen mass loss) are the goals here.

Scatterplots are the most intuitive way to check the dependence structure of a dataset. Figure 6 illustrates the scatterplot matrix for mean area, mean circularity, collagen mass loss, and time for groups CCO and CCT. Other interesting variables (such as the major axis of the fitted ellipse or even mean roundness of the objects in each image) were not included due to their strong dependence on the area or circularity. The area and circularity of the extracellular matrix provide information about cellular growth and shape.

The strongest dependence was found between the mass loss and time variables. In fact, the collagen mass loss strongly depends on time in a relatively complex nonlinear way. It seems that collagen degradation is composed of two main steps that could be related to two main degradation processes. This trend was observed in both groups CCO and CCT, with only slight differences between them.

In addition, an asymptotic type of dependence between circularity and time was observed, whereas mean area seemed to increase with time. This increase seemed to be of the sigmoid type when the CCO group was studied, but in the case of the CCT group, the type of trend was not clear due to high dispersion.

If the dependence structure of the features of the CO group is studied and compared with that of groups CCO and CCT (see Fig. 7A), we note that there are no trends between variables (either linear or nonlinear), and scattering is substantial. This finding is in agreement with the fact that CO is the control group without stem cells.

Furthermore, Fig. 7B depicts the relations between the variables when a logarithmic transform is applied to the mean area and mean circularity. We see that there could be linear relations between the logarithmic transform of the area and circularity as a function of time. In addition, there could be slightly linear relations of the collagen mass loss with the logarithmic transform of the area and circularity. This information will help to estimate the proper statistical models for the degree of collagen degradation.

Collagen degradation modeling

Estimating the underlying model that explains the degree of collagen degradation depending on time is the goal here. This model can provide information about the different steps of degradation and can serve for forecasting tasks. The mass loss feature is the critical variable chosen for characterizing the degree of type I collagen degradation.

According to the results of the exploratory analysis, the use of a sigmoid-type nonlinear function is needed to describe the relation between mass loss and time. Moreover, it seems that the main trend can be the result of the sum of two sigmoid trends.

A nonparametric GAM with penalized regression splines was applied to confirm the results of descriptive analysis. In fact, in Figs. 8A and 8B, the GAM estimates for groups CCO and CCT are shown. Taking into account the fit and bootstrap confidence intervals, we infer that the main trend could be composed of two sigmoid curves, each one corresponding to different degradation processes.

A more accurate fit was obtained for the CCT group and explained 99% of the collagen mass loss information. In any case, it seems that we can obtain accurate estimates of mass loss as a function of time taking into account the high values of the determination coefficient and the relatively narrow confidence intervals.

The next step was to fit a parametric model based on the adequate nonlinear function. The aim was to obtain a transfer function that, in addition to interpretable parameters, could provide predictions and allow us to compare different groups. The sigmoid trend and the fact that the mass loss increases monotonically support the use of logistic functions. In accordance with the GAM results, a nonlinear regression model based on a mixture of two logistic functions (defined by three parameters) is proposed. This type of nonlinear regression model has been successfully applied to separate overlapping degradation processes within the framework of thermal analysis in degradation studies (Tarrío-Saavedra et al., 2014). Figures 8C and 8D shows the logistic-mixture-model fittings that estimate the mass loss depending on time for groups CCO and CCT, respectively. Estimation and prediction asymptotic confidence intervals are provided too. Table 3 shows signification analysis and the model parameters. Signification analysis, fitted trends, confidence intervals, and the achieved goodness of fit (R² = 0.94 and R² = 0.99 for groups CCO and CCT, respectively) justify the use of the logistic mixed model to explain the mass loss of collagen. The first component is characterized by an inflection point at ∼2 days and a saturation asymptote at 50% of collagen mass loss. On the contrary, the second logistic component is defined at an inflection point of ∼16–17 days and an overall mass loss of between the 21–13 of weight percent loss. The meaning of the two logistic components can be guessed taking into consideration control sample CO. In fact, the samples of the CO group lost mass almost immediately, and the mass remained constant at ∼30% mass loss. Thus, this mass loss is related to another factor apart from cellular activity. Deeper research is necessary to identify the nature of the degradation process (either hydrolytic or otherwise). In addition, we infer that the first logistic component depends, to a large extent, on another degradation factor aside from cellular differentiation and growth. In fact, considering the porosity of collagen, the degradation due to cell growth and differentiation should take place mainly after longer periods. Regarding the second logistic component, it may be more related to cell activity (given that the inflection point is ∼16–17 days) that corresponds to ∼20% collagen mass loss after 28 days.

The progression of degradation of type I collagen can be compared between groups CCO and CCT through examination of the logistic components. Figure 8E shows the fitted first logistic component for groups CCO and CCT. There are only slight differences, but it seems that the degradation process tends to begin even before the inflection point μ_CCO(=1.73days) < μ_CCT(=2.47days). The results on bootstrap confidence intervals for parameter µ(implemented using the nlstools package (Baty et al., 2015) with 1,000 resamplings) suggest that CCO samples tend to begin to degrade before the CCT samples do. There are, however, no differences between the groups in terms of the mass loss percentage and the rate of mass loss for this first degradation process. Furthermore, in Fig. 8F, the second logistic components of groups CCO and CCT are compared. The rate of degradation corresponding to the CCO group is significantly greater than that corresponding to the CCT group. This finding can be interpreted via analysis of the scale parameter. In fact, the bootstrap confidence interval for the scale parameter is placed at lower values than those corresponding to group CCT (at lower values of scale, the logistic curve has a higher rate of change). Therefore, the osteogenic medium slightly promotes the degradation of collagen, advancing the degradation process and increasing the rate of degradation.

Cell growth modeling

The mean area and mean circularity of the identified objects of the extracellular matrix obtained from each image at each time point provide important information about cell growth and differentiation. Accordingly, the dependence relation of the area and circularity with time was characterized and statistically modeled. Figure 9 shows the best models fitting the area and circularity variation depending on time, which were chosen according to goodness of fit and parameter signification criteria.

Collagen degradation depending on cell growth modeling

The next step is to determine the dependence relation between collagen degradation (through its mass loss percentage) and cell growth and cell differentiation indicators such as mean area and circularity of the extracellular matrix.

On the basis of the results in the previous sections, logarithmic transformation of the area and circularity was proposed above to obtain more informative, simpler, and more accurate models. Indeed, linear regression models are proposed to explain the collagen mass loss with respect to circularity and the area (after natural logarithms are applied). Figure 10 presents different linear models fitted to the real data obtained from groups CCO (Figs. 10A, 10C) and CCT (Figs. 10B, 10D). Very informative linear models were fitted to estimate the mass loss versus area (massloss =  − 5.217 + 10.578ln(area) with R² = 0.70) and depending on circularity (massloss =  − 19.4 − 14.9ln(circularity) with R² = 0.88) for the CCO group. It is important to stress the strong relation between the mass loss and circularity. The obtained model could be useful even for prediction tasks. The same trends were found in the CCT group, but high dispersion prevented obtaining more informative models. In conclusion, the degree of degradation, measured as the mass loss percent, is proportional to the logarithm of extracellular-matrix area and inversely proportional to the logarithm of cell circularity.

A more complete multivariate linear model that accounts for all the collagen degradation sources of variation was also estimated and validated by a threefold cross-validation procedure. The estimated model expression (see Table 4) is ${\displaystyle \widehat{Massloss}=-17.0396-0.215CCT-2.1751CO-0.0599ln\left(circularity\right)+3.7804ln\left(mode\right)+0.1760ln\left(area\right)}$ where $\widehat{Massloss}$ is the mass loss estimates; CCT and CO are dichotomous variables equal to 1 if the corresponding group is CCT or CO, respectively; and mode is the grayscale histogram mode corresponding to each image. It is characterized by R² = 0.89. The latter variable is included to take into account information about the image color (possibly related to cellular growth) in the model. Increasing the extracellular-matrix area and histogram mode increases the collagen mass loss, whereas decreasing circularity increases the mass loss. The effects in CO and CCT groups are a decrease in the collagen mass loss due to cell activity, when compared with the CCO reference group.

Figure 11A illustrates relative importance of the group, extracellular-matrix histogram mode, circularity, and area for the type I collagen mass loss, in % of R². For this purpose, the R² contribution averaged over orderings of regressors and lmg metrics (Grömping, 2006) were used. Almost all the variability in the collagen mass loss could be explained by group (∼56% of R²) and mean area variation (∼23% of R²) (with preapplied logarithm transformation), but the contributions of circularity (∼8% of R²) and mode (∼14% of R²) are significant too. The overall determination coefficient of the multivariate linear model is 0.89.

Figure 11B shows the results of the cross-validation procedure for measuring predictive accuracy of the multiple linear regression model. The data are randomly assigned to three different folds. During three iterations, each fold is removed (test sample), while the remaining data (training sample) are employed to evaluate the model and then to predict the mass loss corresponding to the removed observations (corresponding to the test group). In fact, observed versus predicted values are depicted in a scatterplot. The cross-validation predictions for each of the three groups in which the data are divided are indicated with different colors. The data points revealed linear trends almost coincident with each other and very close to the bisector. These findings support high predictive accuracy of the newly developed multivariate model for analysis of collagen mass loss. In conclusion, prediction of the degradation degree of collagen (measured as the percentage of mass loss) could be done via the predictors based on image analysis of the extracellular matrix.

Data, workplace and scripts corresponding to the statistical analysis can be found in Dataset S2, Dataset S3 and Code S1.