Analytical model and Monte Carlo simulations of polymer degradation with improved chain cut statistics

The degradation of polymers is described by mathematical models based on bond cleavage statistics including the decreasing probability of chain cuts with decreasing average chain length. We derive equations for the degradation of chains under a random chain cut and a chain end cut mechanism, which are compared to existing models. The results are used to predict the influence of internal molecular parameters. It is shown that both chain cut mechanisms lead to a similar shape of the mass or molecular mass loss curve. A characteristic time is derived, which can be used to extract the maximum length of soluble fragments l of the polymer. We show that the complete description is needed to extract the degradation rate constant k from the molecular mass loss curve and that l can be used to design polymers that lose less mechanical stability before entering the mass loss phase.


Introduction
The design of degradable polymers often relies on cleavable functional groups in the main chain [1][2][3][4]. The availability of water in relevant system environments such as the human body or nature created a special interest in hydrolysable bonds, such as ester, anhydride or orthoester bonds. Macromolecules are split into two smaller fragments in each hydrolysis step. If these polymer fragments are small enough, they become water soluble and the mass of the polymer material decreases depending on the diffusibility of these molecular fragments through the amorphous domains of the partially degraded material. In addition, the diffusion of water in the polymer plays an important role. Bulk erosion occurs homogenously in the entire polymer if water diffusion is fast in comparison to the kinetics of hydrolysis.
Bulk erosion of a homogenous amorphous polymer material starts with an induction phase and continues with a degradation phase. In the first phase, primarily water diffuses into the polymer material. The duration of such a phase depends on device size and geometry. Because water diffusion is much faster than the hydrolysis reaction in case of bulk erosion, there is no significant loss in average chain length or mass of the polymer in this phase. The water concentration in the polymer at the end of the induction phase is equal to its saturation limit, which is about c s ~ 1 wt% for most polymers.
In the following degradation phase, the macromolecules degrade via hydrolysis into smaller fragments. Due to the low concentration of water within the polymer, the concentration of degradable chemical bonds is higher than the concentration of water. That means that every water molecule forms a reaction pair with a degradable bond. If these pairs are inside a reaction volume of radius r, a degradation reaction takes place. The relation between diffusion constant D of water within an amorphous polymer and the reaction constant k defines if the hydrolysis is reaction or diffusion controlled. Hydrolysable polymers with reaction constants of k ~ 10 -2 to 10 -9 s −1 [5,6] and water diffusion constants of D ~ 10 -8 to 10 −10 cm 2 /s [7] have a slow reaction in comparison to the fast diffusion of water within an amorphous polymer, which means that hydrolysis in a polymer bulk is reaction controlled. This can be shown mathematically using the criterion for the transition between a diffusion controlled and a reaction controlled erosion process by Doi et al. [8] This is valid until most of the polymer mass has been degraded into soluble fragments. When these fragments diffuse out of the reaction space, the free volume increases and additional water molecules enter the polymer, which leads to an increase in the water concentration.
As the number of cleavable segments in a chain and the number of chains in a polymer are large, statistical models have to be used to describe reaction controlled polymer degradation. Established models assume that the molecular mass decreases in a single exponential fashion, while a more elaborate version suggests an inverse dependence between molecular mass and time [9,10]. This means that the reaction rate constant can be extracted by means of a linear regression in a suitable plot. However, these models are typically based on the assumption that the number of hydrolysable bonds is constant and are therefore only valid for the very early degradation phase. In reality, the number of breakable bonds is decreasing during the degradation process, so these simplified approaches do not predict the correct shape of the experimental molecular weight loss curve. Here, we hypothesize that a correct quantitative description of the degradation process, which includes the decreasing number of breakable bonds, enables better predictions about the influence of molecular parameters like molecular mass or hydrophilicity on polymer erosion. Therefore, a degradation curve for the molecular mass and for the total polymer mass is derived, that accounts for both the decreasing probability of a chain cut and the chain length-dependent probability for the generation of a soluble fragment. The treatment is limited to amorphous polymers consisting of linear macromolecules with identical repeating units (RU) (homopolymers). A repeating unit is defined as the repeating sequence containing one cleavable bond. Then, the number of cleavable bonds in one chain is the number of RUs minus one. Only the degradation phase is considered, i.e., water is assumed to be present homogenously within the polymer. Furthermore, it is assumed that polymeric fragments with less than l repeating units are solubilized, leading to instantaneous mass loss. Both assumptions are generally valid for small devices such as microparticles. Following these assumptions, the probabilities of chain cleavage for a random chain cut and a preferential chain end cut mechanism are formulated. In both cases, the corresponding system of differential equations (DE) is solved and time dependences for the number of repeating units RU(t) of the partially degraded macromolecules, the number of chain ends NE(t) and the degree of polymerization N P (t) are derived. The analytic calculations are further verified by Monte Carlo (MC) calculations. Next, characteristic quantities of these degradation curves are extracted, which help to interpret specific stages of degradation and from experimentally measured data. Our result is compared to the established first-order approximations for molecular mass degradation and it is shown that these approximations are not able to characterize the decrease of the molecular weight quantitatively correct, which means that the decreasing number of cleavable bonds during the degradation process has to be taken into account for a correct extraction of the degradation rate constant k. Finally, the results are summarized and discussed.

Solution for a random chain cut mechanism
In this section, the degradation of amorphous homopolymers is analyzed. The polymer contains chains with RU(0) repeating units and NE(0) chain ends at the beginning of the degradation.
Firstly, a random chain cut mechanism is considered. In this case, the probability to have a chain cut between two repeating units is independent of the position of these repeating units within the chain. Two chain cuts are possible: A chain cut at the end of a chain at one of the last l degradable bonds on both sides (preferential chain end cut CE), which removes a soluble fragment of length N diff ≤ l from the long chain and a chain cut within the chain farer away than l bonds from any of the two chain ends (inner chain cut IC) that cuts the chain into two smaller chains, which are both longer than l. Here, l is the maximum length of a soluble fragment. A polymer with an average degree of polymerization N P has, on average, NE/2 (N P -1) = (RU -NE/2) cleavable bonds. If N P > 2 l + 1, the probability for a chain cut within the l segments at both ends of the chain is 2 l/(N P -1). In this case, a soluble fragment of length N diff is removed from the chain. On average, a soluble fragment has a length of N diff ~ (l + 1)/2. Similarly, the probability for a chain cut within the segments that are further than l segments away from the chain ends is 1 -2 l/(N P -1). In this case, the number of repeating units is unchanged, but two chain ends are created (see Fig. 1A).
This leads to the following DE for a polymer with N P > 2 l + 1: Here, NE 1 and RU 1 represent the first part for N P > 2 l + 1 of the curves NE(t) and RU(t) with initial conditions NE 1 (0) = NE(0) and RU 1 (0) = RU(0). Similarly, the indices 2 and 3 represent the second part for l + 1 < N P < 2 l + 1 and third part for N P < l + 1 of these curves (see below). Some polymers show a preferred degradation reaction at the end of the chain [11][12][13][14]. A parameter P inner is introduced to represent the relative reactivity of inner chain bonds with reaction constant kP inner and bonds at the end of the chain with reaction constant k. For a random chain cut mechanism, P inner = 1, while for a preferential chain end cut mechanism, P inner < 1.
Independent of the mechanism, if l < N P ≤ 2 l + 1, the chains are so small, that every chain cut creates at least one soluble Article fragment. Chain cuts in the center of the chain can even cleave the chain into two soluble fragments (dissolution cut DC).
The probability of such a chain cut increases with decreasing fragment length and is (2 l + 1-N P )/(N P -1). If such a chain cut occurs, the chain with, on average, 1/2(3 l + 1) repeating units and 2 chain ends is dissolved in water (see Fig. 1B). That means, that the DEs for a polymer with l < N P ≤ 2 l + 1 are Here, NE * (t) = NE 2 (t-t 1 ) and RU * (t) = RU 2 (t-t 1 ) are the time-shifted curves of the second part of the full degradation curves NE(t) and RU(t) with initial conditions NE * (0) = NE 1 (t 1 ) and RU * (0) = RU 1 (t 1 ). N P (t 1 ) = 2 l + 1, i.e., t 1 is the time of transition between the first and second system of DEs. If l < N P ≤ l + 1, every chain cut will create two soluble fragments. That means that all l + 1 repeating units and 2 chain ends are dissolved in water in chain cuts of every of the l degradable bonds. The DE's for l < N P ≤ l + 1 are Here, NE ** (t) = NE 3 (t -t 2 ) and RU ** (t) = RU 3 (t -t 2 ) are the time-shifted curves of the third part of the full degradation curves NE(t) and RU(t) with initial conditions NE * (0) = NE 2 (t 2 ), RU * (0) = RU 2 (t 2 ) and N P (t 2 ) = l + 1. t 2 is the time of transition between the second and third system of DEs.
Single degradation steps in one chain are described by the processes derived above in a stepwise fashion, e. g., the numbers of chain ends and repeating units are changed by a constant number within a chain. Polymers usually contain many chains, which degrade in parallel. The degradation functions NE(t) and RU(t) describe the average number of chain ends and repeating units of all chains and are continuous functions.
The systems of differential equations are solved separately in the Supplementary Information 1-3. The derivation for N P > 2 l + 1 can be found in SI Eqs. 1-25. The solution is for P inner ≥ 4l 2 +4l 4l 2 +4l+1 and for P inner < 4l 2 +4l 4l 2 +4l+1 . The derivation for l + 1 < N P < 2 l + 1 can be found in SI Eqs.26-35. The solution is with with t 1 as the time when N P (t 1 ) = 2 l + 1: for P inner ≥ 4l 2 + 4l for P inner < 4l 2 + 4l This is the solution for ω ≥ 0. The solution for ω < 0 is: The derivation for l < N P < l + 1 can be found in SI Eqs. 36-44. The solution is The full functions NE(t), RU(t) and N P (t) = RU(t)/(NE(t)/2) are plotted in Fig. 2 for the degradation of one chain with initially RU(0) = 100 repeating units and NE(0) = 2 chain ends and typical values for the maximum length of soluble fragments l = 4 and reaction rate constants k = 10 -3 . A value of l = 4 is typical for PCL. l is determined experimentally as the highest molecular weight fragment found in the aqueous phase. The reaction rate constant is typically measured in reciprocal time units, i.e., per second, per hour, per day etc. The unit of time for such a degradation experiment would be the corresponding timescale, i.e., seconds, hours or days. For example, a reaction rate of k = 10 -3 h −1 can be regarded as typical for poly[(rac-lactide)-coglycolide] (50 : 50) in water at 37 °C and pH = 7 [15]. Accordingly, the x-axis would be in hours. The model predicts 50% mass loss after 350 h, which is in agreement with the experimental observation [15].
The number of chain ends increases from its initial value to its maximum value NE(t 1 ). The degree of polymerization decreases, and the slope of the decrease is determined by the increase of NE(t) from chain cuts occurring too far away from the chain end to create soluble fragments. Polymer mass loss starts at the moment when the chain length starts to decrease, but the trends are opposing: While molecular mass loss is fastest at the beginning, mass loss is initially extremely slow. This means that in degradation experiments, a constant molecular mass can only be observed in the induction phase, where water enters the polymer, which is not considered here. Polymer mass loss is accelerated for t ≥ t 1 . The degree of polymerization decreases slowly until l + 1 at t 2 . Behind t 2 , the low number of remaining hydrolysable bonds leads to a slow exponential decrease of NE(t) and RU(t) and N P varies around l + 1 until the last polymer fragment is dissolved. When neglecting the increasing probability to generate water soluble fragments for t > t 1 , the full degradation curves are still quite well replicated (dotted violet line in Fig. 2), but with a slower decrease of the polymer mass. Figure 3 shows the degradation curves for different initial number of chain ends NE(0) (Fig. 3A), initial number of repeating units RU(0) (Fig. 3B), lengths of soluble fragments l (Fig. 3C) and hydrolysis constants k (Fig. 3D) for a random chain cut mechanism. With increasing initial number of chain ends NE(0), but constant polymer mass, the degradation curves of NE(t) and RU(t) are shifted to shorter times and the maximum number of chain ends increases slightly. This situation is similar to the degradation of one long chain, which has been already degraded by some inner chain cuts, but no chain end cuts. The maximum number of chain ends also increases with increasing NE(0). Interestingly, the maximum number of chain ends only increases insignificantly from NE(t 1 ) = 16.40 for NE(0) = 2 to NE(t 1 ) = 17.76 for NE(0) = 14. This shows that the molecular mass loss is not much affected by the initial molecular mass.
If the initial polymer mass is increased, the maximum number of chain ends increases significantly. This increase is nearly linear. For example, the maximum number of chain ends increases from NE(t 1 ) = 37.03 for RU(0) = 50 to NE(t 1 ) = 368.65 for RU(0) = 500. Interestingly, the time to reach this maximum is nearly unchanged. The time to degrade to N P (t 1 ) = 2 l + 1 is nearly independent of N P (0). The time t 1 can therefore be used as a characteristic property of the degradation curves (see below).
Two important parameters, which characterize different degradable polymers, are the hydrophilicity of their monomers and their hydrolytic stability. The former is related to the maximum length of their soluble fragments l, while the latter is indicated by their degradation rate k. A higher value of l causes a faster polymer mass loss because degradable fragments are produced with less chain cuts. Similarly, a higher hydrolysis rate k, e.g., due to catalysis or higher water concentration, leads to a faster mass loss of the polymer. However, the influence of k and l on the molecular mass is opposite: An increase in l leads to a slower decrease of the molecular mass, while an increase in k leads to a faster decrease of the molecular mass. This can be explained with the number of chain ends NE(t) or the total number of chains during the degradation. While a change in k only changes the time scale of the degradation functions NE(t) and RU(t), but not its absolute values, a change in l is related to a change in the maximum number of chains at t 1 . For example, NE(t 1 ) decreases from 29.64 for l = 2 to 7.05 for l = 14, so the dependence is almost linear. As molecular mass loss is the primary cause for reduction of elastic moduli or tensile strength, increasing l allows to design polymers that are mechanically more stable when they reach the phase of mass loss. Since l is a threshold for solubility, it relates to hydrophilicity. A high value of l is, e.g., observed in polymers with a high ratio of polar to unpolar groups like PLGA.

Solution for a preferential chain end cut mechanism
If the degradation happens preferentially at the end of a chain, a chain end cut mechanism with P inner < 1 for a cut in the inner part of the chain is used to describe the degradation. The solution for a chain end cut mechanism is plotted for different values of P inner in Fig. 4A. Lower values of P inner lead to a slightly slower increase of NE until t 1 . In general, the degradation curves are quite unaffected by a lower reactivity of inner chain cuts. This means that the observed degradation rate constants correspond to the degradation constant of the terminal bonds of a chain and that a chain end cut mechanism instead of a random chain cut does not change characteristic degradation parameters (see below).
Monte Carlo simulations for N polymer = 95 independently degrading polymers have been performed. Figure 4B shows that the averaged number of chain ends and repeating units and the degree of polymerization are similar in shape and amplitude to the analytical solution, which proves that the analytical solution describes indeed the statistical degradation of a polymer with a random and a preferential chain end cut mechanism. A decrease of N P below l is not possible because every degraded polymer leads to an increase in the average degree of polymerization of the remaining ensemble.

Characteristic Quantities
In the previous sections, the degradation curve for a reaction controlled degradation is discussed. It is shown that a chain end cut mechanism does not change the general shape of the degradation curve over a random cut mechanism. In this section, characteristic quantities are extracted from the degradation curve of a random chain cut mechanism, which are helpful to characterize and interpret experimentally measured degradation curves. As discussed above, the time t 1 , at which the number of chain ends reaches its maximum, is nearly independent of the initial chain length or the number of molecules. The number of repeating units for large initial chains (RU(0) > > NE(0)) at t 1 is: Figure 5 shows the value of this limit as a function of l. The expression reaches a constant plateau value for increasing l at ~ 0.736, which is already reached for l = 3 with only slightly higher values for l = 1 (0.75) and l = 2 (0.741). That means that the degradation curve reaches t 1 when ~ 26.4% of the initial polymer mass is degraded. At this point, the degree of polymerization reaches the value of N P (t 1 ) = 2 l + 1. That means that. 1. the measured polymer mass can be used to extract t 1 , the time when the experimentally not directly accessible number of chain ends reaches its maximum and the mass loss is accelerated, and 2. the measured molecular mass at t 1 can be used to extract the value of l for a specific polymer.
This is especially important for block (co)polymers with monomeric units of different l. If the values of l are known for all monomeric units, the value of l of a mixture can be used to extract the relative contribution of the individual blocks.

Comparison to established description
First-order differential equations of the molecular mass lead to solutions, which are either single exponential [16] or linear in 1/M N [9,10] . It is common practice to extract the reaction rate constant from a linear regression of the logarithm or the inverse of the molecular mass as a function of time. Figure 6 shows the logarithm and the inverse of the degree of polymerization from the Monte Carlo simulation with P inner = 1. Both functions are nonlinear. This means that the hydrolysis rate k cannot be directly extracted via linear regression. The result will always vary strongly depending on the position of the degradation curve where the regression is carried out. The degree of polymerization N P (t) = 2RU(t)/NE(t) contains hyperbolic sine and cosine functions in the nominator and in the denominator, which both contribute to the decrease of the molecular mass curve. The leading terms of the Taylor expansion of this function can only replicate a very small part of the curve. It is therefore necessary to fit the decreasing part of the degree of polymerization (or molecular mass) curve with the full solution of the first phase of degradation to extract the correct hydrolysis reaction constant k.

Conclusions
Differential equations for the reaction controlled bulk degradation of amorphous homopolymers have been derived and solved. These equations are valid for the degradation phase following an initial phase required for water uptake, of which the duration depends on device size and geometry. The differential equations are based on chain cleavage statistics of linear chains and consider the decreasing reaction rate caused by decreasing chain length. The solution yields a sigmoid shaped curve for polymer mass loss. The degradation curves are discussed for different values of the initial molecular length, degree of polymerization of the polymer, fragment solubility limit and for different hydrolysis rates. Three regions based on the degree of polymerization are considered separately. It is shown that the transition between these regions can be Figure 6: Degree of polymerization with logarithmic axis (top) and inverse (bottom) of the degree of polymerization from the Monte Carlo simulation with P inner = 1 after an induction phase of t i = 50, which represents the phase, in which water enters the polymer before degradation starts. used to extract important parameters like the length of soluble fragments l from experimentally measured polymer mass loss and molecular mass loss curves. The solution is compared to approaches used in the literature, which use firstorder approximations of the molecular mass to describe the degradation of polymers. It is shown that plotting the logarithm or the inverse of the molecular mass as a function of time is not appropriate to extract the hydrolysis rate k. It is shown instead, that the correct mathematical description of the degradation process is more complex and requires the full solution to extract rate constants from experimental degradation curves. This complete description can be used to predict the degradation of polymers with known k and l and for the design of new materials.

Materials and methods
Monte Carlo (MC) simulations are performed to describe the degradation of polymer materials. Polymers are considered as linear chains of initially RU(0) repeating units and NE(0) chain ends. At every MC step, two RUs are separated with probability p = RU(t)-NE(t)/2. Random chain cuts, chain end cuts or dissolution cuts are performed with the probabilities and step sizes described above. The MC simulation is stopped when all cleavable bonds have been degraded. The simulation is performed 95 times for P inner = 1 for a random chain cut mechanism and P inner = 0.7 for a chain end cut mechanism. Lower values of P inner have not been used because it has been found that even in case of significant chain end cuts, degradation in the main chain occurs in parallel [12].

Funding
Open Access funding enabled and organized by Projekt DEAL. This work was financially supported by the Helmholtz Association of German Research Centers through programoriented funding.

Availability of data and material
Our manuscript has data included as electronic supplementary material, additional data will be made available on reasonable request.

Conflict of interest
The authors declare no competing financial interest.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.