Degradation modeling of degradable copolymers for biomimetic scaffolds

Biomimetic scaffolds provide a suitable growth environment for tissue engineering and demonstrate good potential for application in biomedical fields. Different-sized copolymerized biomimetic scaffolds degrade differently, and the degradation rate is affected by the copolymerization ratio. The study of the degradation property is the foundational research necessary for realizing individualized biomimetic scaffold design. The degradation performance of polyesters with different copolymerization ratios has been widely reported; however, the modeling of this performance has been rarely reported. In this research, the degradation of copolymers was studied with multi-scale modeling, in which the copolymers were dispersed in a cellular manner, the chain break time was simulated, and the chain selection was based on the Monte Carlo (MC) algorithm. The probability model of the copolymer’s chain break position was established as a “roulette” model, whose probability values were estimated by the calculation of the potential energy difference at different chain break positions by molecular dynamics that determined the position of chain shear, thereby fully realizing the simulation of the chain micro-break process. The diffusion of the oligomers was then calculated using the macro diffusion equation, and the degradation process of the copolymer was simulated by three-scale coupling calculations. The calculation results were in good agreement with the experimental data, demonstrating the effectiveness of the proposed method.


Introduction
With the development of tissue engineering in recent years, biomimetic scaffolds, aside from playing a supporting role, can now provide a beneficial environment for cell nutrition, growth, and metabolism, and can provide the technical means for tissue regeneration. Their gradual degradation and absorption of biomaterials decrease the effectiveness of the mechanical properties of the biomimetic scaffolds. At the same time, implanted cells continue to proliferate and secrete extracellular matrices in vivo, eventually forming the corresponding tissues or organs to repair and achieve reconstruction after trauma. The rate of decreasing effectiveness of the mechanical properties of the biomimetic scaffolds should be consistent with the rate of increase in tissue regeneration. However, the autocatalytic hydrolysis reaction leads to a size dependency [1]: the degradation rate in the center of the device is higher than that around the boundaries, which means that larger scaffolds may degrade faster than smaller scaffolds. Therefore, the degradation theory Copolymer (such as poly(lactic-co-glycolic acid) (PLGA)) biomimetic scaffolds are becoming increasingly widely used in medicine [2−4] because the rates of change of the mechanical and degradation properties can be controlled. Different copolymerization ratios correspond to different properties [5−8] and change during the degradation process, so it is important to specifically study the degradation process of copolymers with different copolymerization ratios. The degradation process of copolymers is complex, with macroscopic changes in molecular weight, mass, and mechanical properties caused by microscopic chain scission (through the hydrolysis reaction), autocatalysis of the degradation process, recrystallization, and diffusion of oligomers. The method of experimental study and requires long periods of time, large amounts of manpower, material, and financial resources, and the size dependency caused by the autocatalysis of the degradation process makes the degradation rate of different equipment sizes unpredictable. Computer modeling is an appropriate supplementary research method that can simulate the degradation process and provide quantitative data, which provides a basis for the optimization of copolymer equipment. Modeling studies of copolymer degradation have not been reported in public literature as most of the current research focuses on modeling the degradation of homopolymers. The various scales of model processing can be divided into single-scale macroscale modeling [9−11] (in the commonly used deterministic method for continuous modeling), microscale modeling [12−14] (used in the common stochastic method of modeling), mesoscale modeling [15,16] (for the commonly used discretization modeling), two-scale models [17], and multi-scale models [18]. However, most of these models are used in the research of the degradation process of homopolymers, and they do not model the copolymer components in the degradation process, nor do they involve the scission study of different copolymer chains. In this paper, the chain scission model of copolymers is proposed. The scission process of the copolymer chain is recorded by the position of the chain scission, and the multi-scale model coupling is carried out to simulate the degradation process of the copolymer, which provides new methods and concepts for the design and performance optimization of copolymers.

Chain shear modeling of copolymers
Microchain scission modeling is often simulated by the dynamic Monte Carlo (MC) method [17,18], but this only determines the time and the specific chain that undergoes the chemical reaction of chain scission (the MC method randomly selects a certain molecular chain at a certain time to break), while the specific location of the break in the chain is not considered. However, the position of the chain scission (end scission or random scission) has a significant effect on the degradation process. Only the conditions under which the end scission and random scission occur are often assumed in these simulations, but this paper proposes a probabilistic model of the chain scission position of the copolymer. The method of calculating the potential energy difference by molecular dynamics provides a theoretical basis for identifying the possible position of the chain scission, and thereby simulates the microscopic scission of the chain more accurately. It is assumed that the molecular chain of the copolymer is continuously distributed as shown in Fig. 1. We use the Materials Studio software to construct the PLGA chain, as shown in Fig. 2.
The potential energy E 0 of the chain is calculated according to the principle of energy minimization. Assuming that the chain breaks at a specific position, for example, a fracture at q 1 , we then recalculate the potential energy of the chain after fracture, denoted as E q1 ; by computation, the potential energy for different chain breaks is calculated, denoted as E qi (i = 1, 2,..., M). The potential energy difference between cases of chain fracture at different positions can be calculated through Eq. (1): Taking the potential energy difference as a theoretical reference, the probability P qi of chain fracture location is found through Eq. (2).
In this way, the position probability model for chain shear fracture is constructed. The probability obtained by calculating the potential energy difference is recorded by means of a probability roulette, and the size of the roulette area is determined by the probability value of P qi , shown in Fig. 3.
Starting from the q 1 position, the anticlockwise direction represents the fracture probability of q 1 , q 2 , …, q M positions in the chain shown in Fig. 1. The roulette programming algorithm is used to select the position of the chain scission, q j (j ∈ [1, M]), which is recorded. The molecular chain breaks from position q j , and the two new chains generated overwrite the position of the bond.

A simplified algorithm for determining the location of chain breaks during degradation
The potential energy difference calculation can be used as the theoretical basis for the chain scission probability. However, in the real degradation process, the chain breaking reaction is so fast that it is unfeasible to calculate the potential difference for each new chain update. In this paper, using the potential energy difference calculation as the theoretical basis, polynomial linear regression is used to construct the chain scission

Copolymer degradation process calculation method
The micro chain scission model is combined with the use of mesoscopic cells and the macro diffusion equation. The degrading polymer can be considered as a heterogeneous material at the mesoscopic scale made up of amorphous, crystalline, and hollow phases. The material is divided into β × β grids referred to as cells. The chain scission reaction happens in each cell at a mesoscopic scale. Cellular states change during degradation, and the molecular weight, chain number, and the proportion of copolymerization components are recorded in the individual cells. The oligomers produced diffuse out of the polymer matrix according to Fick's second law. The molecular weight distribution, chain number, proportion of copolymerization components, and weight loss that vary with time are obtained by the statistics of cellular record. The specific steps are as follows: Step 1: The copolymer PLGA is dispersed into β × β equal cells that are then initialized, the cell (a, b) denoted as Cell[a] [b]. The molecular chain at time t is denoted as Chain_c(c = 0, 1,2,...,H, where H is the total number of chains), containing X PLA ester bond units and Y PGA ester bond units. Step The relationship between neighbor cells is based upon the von Neumann 4-neighborhood model.
Step 3: Traversing all cells, the number of chains C I (I = 0, 1,2,...,M is the number of ester bonds and M is the largest number of ester bonds) is counted.
Step 4: The time interval Δt 1 (designated as the chain scission time step) and specific type of reaction u of the hydrolysis chain scission reaction are determined according to the MC method. In the MC scheme, the ester bond μ is selected for chain breakage according to with the following length of time integration step: in which v = 1, 2, …, M, and X pv , X ol , and X W are the reactant molecules (polymer with the ester bond unit v, oligomer, and water molecule, respectively) [17]. π v1 and π v2 are the rate constants for the non-catalytic and autocatalytic hydrolysis reactions of ester bonds v, respectively. Although the description of this equation considers the molecule degree, the numerical values of the dimensionless rate constants are the same as the reactions constants k 1 and k 2 [17].
Step 5: According to Eq. (2), the fracture position q i of the chain is selected as the location where the chain breaks, and the new number of chains and the positional arrangement of the new chain are then recorded.
Step 6: Accounting for the products of the above microscopic reactions and mesoscopic cells, the physical properties (molecular weight, chain number, weight loss, proportion of copolymerization components, etc.) of the cells are recorded and the cell states are updated.
Step 7: The time value is updated to t 1 = t 1 +Δt 1 . If t 1 < Δt 2 (where Δt 2 is the diffusion time step), we return to Step 3, otherwise we proceed to Step 8.
Step  (7) in which C e is the concentration of the copolymer ester bond, C ol is the concentration of the oligomers, k 1 is the rate constant of the hydrolysis reaction without catalysis, k 2 is the rate constant of hydrolysis reaction with catalysis, and D is the diffusion coefficient described by the equation from Han and Pan [17]: where D 0 is the diffusion coefficient of the oligomer in the polymer, D 1 is the diffusion coefficient of the oligomer in the pores, and ε is the porosity, given by is the total number of copolymer cells.
Step 9: The results of the calculation are produced, recording molecular weight, mass loss, and proportion of copolymerization components, etc. The number of diffusing oligomer is deduced in preparation for the next chain scission reaction.

Calculation
Calculation of the copolymer degradation was based on experimental data with two different copolymerization ratios (PLGA (53:47) and PLGA (75:25)) [19]. To prepare the experimental sample, the polymer was placed in a melt indexer to be melted completely over the necessary equilibration time, then extruded into a rod with a 2.0 mm diameter mold by a piston loaded at high temperature (140 °C ). The degradation experiments were carried out in a phosphate-buffered solution (pH = 7.4) at 37 °C , and the detailed procedures of experiment and measurement are described in Ref. [19].
The parameters of the calculation model were set as follows: the cell dimension β = 1000, the chain scission probability was calculated by taking PLGA (53:47) as the sample to calculate the potential energy difference, and a PLGA copolymer consisting of 113 lactic acid (LA) monomers and 101 glycolic acid (GA) monomers was constructed. The calculation results are shown in Fig. 5. It can be seen that the potential energy difference between the two ends of the chain is larger than the sections in the chain between them, which means that the two ends of the chain are easier to break. The calculated potential energy difference is inputted into the chain scission roulette probability model to screen for the breaking position of the chain. This position, combined with identification of the chain that precipitated the reaction at the time determined by the MC model when the chain break reaction occurred, models the complete detailed simulation process of the microscopic chain fracture reaction. Combining the calculation of microscopic chain scission with the mesoscopic cell approach and considering the diffusion of macroscopic oligomers, the model parameters are shown in Table 1.
The calculated values of the PLGA (53:47) and PLGA (75:25) models are compared with the experimental data as shown in Figs. 6 and 7, respectively. Table 1 Parameters of different copolymerization ratio calculation model.

Autocatalytic treatment
The calculation of the molecular chain scission position described above is based on the potential energy difference between pre-breaking and the post-breaking states. However, the real degradation process is complicated, especially when the oligomer gathers in the matrix to form an acidic environment that results in autocatalytic degradation, in which determination of the chain scission position also becomes more complicated. According to the theory behind the potential energy difference, this model uses the evolutionary rules of cellular automata to account for the autocatalytic phenomenon, that is, when the oligomers occur in a particular cell or its neighbor cells the probability of the cell's random fracture increases, which modifies the roulette model. When the autocatalytic reaction occurs, the probability of random chain fracture increases, while the molecular weight begins to decrease significantly.

Variations in copolymerization composition
The different possible positions of chain scission during degradation lead to a variation in the overall proportions of GA and LA, and variation in the composition of the entire complex. With our model, we calculate the proportion of GA in the above example (LA:GA = 75:25) as a function of degradation time.
The calculation results are shown in Fig. 8, which shows that the proportion of GA does not vary significantly with the degradation process, although the molecular weight decreases significantly (see Fig. 7(a)). The oligomer aggregation causes the autocatalytic reaction to accelerate the degradation until the reaction time reaches approximately 40 days, and results in an evident decrease in the proportion of GA and a significantly increased mass loss (see Fig. 7(b)). The calculation results fit the experimental results well, so we conclude that the calculation model can provide the quantitative data needed to show the degradation process more clearly.

Model sensitivity analysis
The sensitivity of the model should be analyzed with respect to its theoretical predictions. The root mean square error (RMSE) for the normalized average molecular weight and the computation time with the change of model parameter β are used to analyze the model's sensitivity, which is shown in Fig. 9. It can be seen that the model is convergent, and the size of cellular grid, β, affects the computation time. For larger β values, the program requires more memory for the larger cells, and selecting the cells and chains Fig. 9 The sensitivity analysis of root mean square error (RMSE) for normalized average molecular weight (a) and computation time (b) with the change of model parameter β, the experimental data is reproduced with permission from Ref. [19], © Elsevier 2009.
to participate in the chain scission reaction will also require more time. However, larger values of β will yield a more accurate modeling result. Based on this analysis, it can be seen that the choice of parameter β is based on the compromise between time cost and required accuracy. In this paper, compared to the experiment in Ref. [19], β is set as 1000.

Conclusions
In this paper, the degradation process of copolymers is studied by a multi-scale micro chain scission model using mesoscopic cells and the macro diffusion equation to model chain breakage and oligomers diffusion to obtain molecular weight, chain numbers, proportions of copolymerization components, and weight loss. The chain break probability model is used to determine the break position of the chain and uses the MC method to construct the micro chain scission model. The roulette model probability value is calculated by molecular dynamics on the basis of the potential energy difference and is combined with the evolutionary rules of cellular automata for autocatalytic chain break simulation. The microscopic chain fracture model is more effective for simulating chain scission than using the MC method alone; the calculation results of the model are compared with the experimental data of the example copolymer and found to be in agreement, which indicates the validity and effectiveness of this modeling method. The modeling of copolymer degradation in this manner provides a new concept for further research, new methods, and reference for the design and performance optimization of copolymers.