Cellular Automaton Simulation for Degradation of ... - ACS Publications

Subscriber access provided by UNIV OF CAMBRIDGE

Modeling and Informatics Tools

Cellular automaton simulation for degradation of poly lactic acid with acceleratable reaction-diffusion model Chao Guo, and Yi Niu ACS Biomater. Sci. Eng., Just Accepted Manuscript • DOI: 10.1021/acsbiomaterials.9b00015 • Publication Date (Web): 06 Mar 2019 Downloaded from http://pubs.acs.org on March 18, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Biomaterials Science & Engineering

Cellular automaton simulation for degradation of poly lactic acid with acceleratable reaction-diffusion model Chao Guo*, Yi Niu

School of Materials Science and Engineering and Jiangsu Key Laboratory of Advanced Metallic Materials, Southeast University, Nanjing, 211189, China

KEYWORDS: PLA; degradation; cellular automaton; acceleratable reaction; diffusion

ABSTRACT: Biodegradability is a fundamental property of poly lactic acid (PLA). Numerous simulation algorithms based on reaction-diffusion model have been utilized to analyze the degradation behaviors of PLA. In the present work, a cellular automaton (CA) algorithm combined with an acceleratable reaction-diffusion model and coarse-grained kinetic Monte Carlo method is applied to simulate PLA’s degradation behaviors, such as random or end scission and crystallization of PLA chains, diffusion of soluble oligomers.

ACS Paragon Plus Environment

1

ACS Biomaterials Science & Engineering 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 42

The CA algorithm can reveal the global changes of molecular weight, mass, reaction number and soluble oligomers’ number generated by hydrolysis as well as the local distributions of molecular weight, soluble oligomers’ number and cellular state. The calculation result of experiment demonstrates that such a CA model can accurately simulate the change of molecular weight and mass loss simultaneously. The effects of hydrolysis mode, reaction rate constant, diffusion coefficient, device size and pore structure on PLA’s degradation behaviors especially the change in the molecular weight and the famous autocatalysis effect are comprehensively investigated. A novel classification method of molecular weight change curve is presented and the effects of various factors are concluded. In general, big reaction rate constant, small diffusion coefficient, big device size, solid or low porosity and none or few initial soluble oligomers can more easily generate type I curve which corresponds to a strong autocatalysis hydrolysis.

INTRODUCTION


2

Page 3 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Poly lactic acid (PLA) is a famous polyester biomaterial and widely utilized as surgical suture, bone fixing material, drug delivery carrier, tissue engineering scaffold, and so on 1-4.

The biodegradability of PLA is the key of these applications, therefore, the research

of PLA’s degradation behaviors is important and popular 4-9. The most frequent method to investigate the degradation of PLA is immersion experiment in vitro with some special solution such as phosphate buffered solution (PBS) or simulated body fluid (SBF). Numerous factors involving PLA property (such as its molecular weight and distribution), device property (for example, the size and porosity of PLA device) and surrounding condition (including surrounding temperature, solution type, static or flowing solution, etc.) were discussed by different researches

10-16.

However,

some disadvantages limit the effect of degradation experiment. First, experimental results themselves only reveal superficial conclusion and the in-depth analysis of these results needs the guidance of degradation theory or mathematical model. Second, the longer degradation period and massive impact factors make it difficult to fully understand the degradation behaviors of PLA by experiment during a short time.


3


Page 4 of 42

Because of its quick calculation speed and the ability of combining with various kinetic equations, computer simulation has been developed to discover PLA’s degradation in recent years. Several approaches including molecular dynamic (MD), Monte Carlo (MC), cellular automaton (CA) and finite element method have been presented to simulate the degradation process of PLA

17-24.

In these simulations, reaction model based on

autocatalytic and non-catalytic reaction dynamic equations is the most common mathematical model and can obtain the overall evaluation of PLA, such as its molecular weight change and degree of crystallinity change. Diffusion model according to Fick’s law has been frequently added to the simulation to explain the special autocatalytic effect of PLA 19-20, 25-32. Other simulation methods have also been joined into the reaction-diffusion model. Kinetic Monte Carlo (KMC) simulations with random scission and end scission have been applied by several research groups to evaluate PLA’s mass loss as well as the drop of PLA’s molecular weight

17, 20, 23.

CA approaches have been employed in

different literatures to explore the local degradation details 18, 23. Recently, erosion models have been utilized to match the mass loss between the simulations and the degradation


4

Page 5 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


experiments

23-24,

suggesting that the traditional reaction-diffusion model is not

appropriate to simulate the mass change of PLA. Analyzing the effects of various reaction mechanisms and corresponding parameters on PLA’s degradation behaviors is an interesting work and can make a deep insight of PLA’s degradation. Han et al. 20 simply estimated the effect of PLA chain cleavage mode on the production of soluble oligomer during PLA’s degradation process. Gleadall et al. 30-31

gave an analysis scheme to reveal the effects of hydrolysis mechanism, PLA chain

cleavage mode, reaction rate constant, initial molecular weight and residual monomer, but the diffusion of soluble oligomers was ignored in part of the analysis process. In the present paper, a comprehensive CA simulation combined with an acceleratable reaction-diffusion model and KMC method is presented to match the mass loss of PLA during the degradation as well as the change in PLA’s molecular weight. The effects of different reaction mechanisms and numerous additional factors promoting or preventing autocatalysis effect on PLA’s degradation behaviors are systematically analyzed in the simulation work. Local degradation details of various conditions are also investigated.


5


Page 6 of 42

THEORETICAL BASIS Simulation hypotheses The necessary hypotheses of our CA model are as following: 1. PLA is degraded by bulk erosion since surface erosion only happens in special condition such as very high surrounding pH value 5, 33. 2. The degradation is initially taken place in the amorphous zone of PLA. Crystallinity zone will not degrade before the degradation of the amorphous zone is finished 34. 3. Chain cleavage of PLA can be generated by random chain scission reaction or end chain scission reaction. New insoluble PLA chain yielded from old chain cleavage may crystallize with a probability of p (between zero and one)

23, 35

or changes into a soluble

oligomer when its ester bond number is less than seven 36-37. 4. The diffusion rate of water is faster than the reaction rate of PLA chain cleavage, so water is sufficient in the degradation 20, 25. 5. The space of surrounding is big enough and the diffusion rate of soluble oligomer in water is faster than oligomer’s diffusion rate in PLA, so the concentration of soluble oligomer in the surrounding is set to zero during the whole degradation process.


6

Page 7 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. A special number average molecular weight threshold Mn,th exists in the degradation, when the number average molecular weight Mn of a zone is under Mn,th, an accelerated hydrolysis reaction happens in this zone. Cellular states In our CA model, each cellular has one of the following two states: pore or material. At the initial stage, only the surrounding or pore structure of porous PLA device has a pore state while all PLA are assigned as the material state. During degradation, the state of a cellular is changed from material to pore when all PLA chains change to soluble oligomers. Degradation dynamic equations and diffusion equations The PLA chain cleavage reaction obeys the following dynamic equations 𝑑𝑅random/𝑑𝑡 = 𝑘ac𝑘rn𝐶a ― ester + 𝑘ac𝑘ra𝐶a ― ester[𝐶oligomer/(1 ― 𝑋𝐶)]1/2

(1)

𝑑𝑅end/𝑑𝑡 = 𝑘ac𝑘en𝐶a ― end + 𝑘ac𝑘ea𝐶a ― end[𝐶oligomer/(1 ― 𝑋𝐶)]1/2

(2)

in which Rrandom and Rend represent mole concentrations of random scission and end scission, krn, kra, ken and kea are the kinetic rate constants for non-catalytic random scission, autocatalytic random scission, non-catalytic end scission and autocatalytic end scission, respectively, Ca-ester, Ca-end and Coligomer represent mole concentrations of amorphous ester bonds, amorphous end chains and soluble oligomer’ ester bonds, XC is the degree of crystallinity, while kac is acceleration coefficient and can be represented as


7


𝑘ac =

{

Page 8 of 42

1, 𝑀𝑛 > 𝑀𝑛,𝑡ℎ 𝑀 𝑎exp ( 𝑛,𝑡ℎ 𝑀𝑛 ― 1), 𝑀𝑛 < 𝑀𝑛,𝑡ℎ

(3)

where Mn is number average molecular weight of a material cellular, a is a constant. The model can turn into the traditional reaction-diffusion model if threshold Mn,th is set to zero or change into the erosion model if a is equal to infinity. According to Han’s coarse-grained KMC simulation algorithm

20,

the corresponding random

scission reaction number Nrandom and end scission reaction number Nend in a time interval Δt can be approximately expressed as 𝑁random = 𝑑𝑅randomΔ𝑡/𝑑𝑡

(4)

𝑁end = 𝑑𝑅endΔ𝑡/𝑑𝑡

(5)

The diffusion equation of soluble oligomer can be expressed by Fick’s second law 𝑑𝐶oligomer/𝑑𝑡 = 𝑑𝑅oligomer/𝑑𝑡 + div[𝐷grad(𝐶oligomer)]

(6)

in which Roligomer represent mole concentration of soluble oligomer, D is the diffusion coefficient of oligomers and can be calculated by 38 𝐷 = 𝐷PLA + (1.3𝜀2 ― 0.3𝜀3)(𝐷pore ― 𝐷PLA)

(7)

where DPLA and Dpore are diffusion coefficients of oligomers in the amorphous PLA and liquidfilled pore, while ε is the porosity of the cellular and can be determined by (8)

𝜀 = 1 ― 𝐶ester/𝐶ester,0

in which Cester and Cester,0 represent current mole concentration and initial mole concentration of total ester bonds including amorphous and crystallized PLA chains in the cellular, respectively. Simulation process


8

Page 9 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The whole simulation process based on above hypotheses and equations can be divided into the following steps. Step 1. A two-dimensional CA model with a size of M×N is established and the state of each cellular was defined by its feature (surrounding, material and pore). Step 2. Corresponding cellular properties including the number and distribution of PLA chains with various molecular weights are assigned for each material cellular depending on the logarithmic normal distribution of weight average molecular weight

18.

Step 3. Nrandom and Nend in the time interval Δt are determined according to equation (4) and equation (5), the corresponding random chain scissions and end chain scissions are simulated one by one and the average molecular weight are calculated for each material cellular at time t. Especially, if the number average molecular weight of a material cellular is lower than Mn,th, its acceleration coefficient no longer maintains one and calculates using equation (3), that is, the cellular undergoes an accelerated degradation. Meanwhile, a random number r between zero and one is generated for every new insoluble PLA chain after chain scission. If r is less than default crystallization probability p, the chain changes to a crystallized chain. In addition, oligomers’ mole concentrations Roligomer before and after


9


Page 10 of 42

all PLA chains’ scissions are gained by counting the number of soluble oligomers and then dividing the corresponding volume to calculate the difference value of dRoligomer /dt, i.e., ΔRoligomer/Δt. Step 4. A global oligomer diffusion is simulated using alternating direction implicit finite difference method for the entire CA model according to equation (6). After diffusion, local CA parameters including cellular average molecular weight, mass, PLA chain number and oligomer number and global CA parameters such as the number of material cellular, average molecular weight, total mass, PLA chain number, oligomer number and the degree of crystallinity are counted, respectively. Total mass m is calculated by counting the mass of every PLA chain remaining in the model according to its molecular weight and adding them while degree of crystallinity XC is determined by counting the molecular weights of crystallized PLA chains and total chains and dividing the latter by the former. Then, mass loss is gained by subtracting current total mass m with original total mass m0. A schematic diagram of PLA chain scission reaction and oligomer diffusion process happened in Step 3 and Step 4 is given in Figure 1.


10

Page 11 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Step 5. Degradation time increased Δt. The simulation is finished if the degradation time reaches the maximum simulation time or all PLA chains are degraded into soluble oligomers. Otherwise, the simulation returns to Step 3.

Figure 1. Schematic diagram of ordinary or accelerated reaction and diffusion process from time t to time t + Δt RESULTS CA model validation


11


The experimental results about size-dependence by Grizzi et al.

Page 12 of 42

10

are simulated and

shown in Figure 2 to affirm the validation of above CA model. Primary simulation parameters are listed in Table 1. Zhang et al. 23 and Sevim et al. 24 have simulated the result of the 2mm PLA plate using different erosion models, Zhang’s model is more consistent with the result of molecular weight change, while Sevim’s simulation matches mass loss well. Compared with these erosion models, our simulation accuracy is similar to that of Zhang’s molecular weight change and Sevim’s mass loss. Moreover, our mass loss result about the 0.3mm PLA film is worse than Sevim’s simulation while the change in molecular weight is more effectively simulated. In other words, our simulation is on the same level with Sevim’s calculation result. All these facts indicate that our acceleratable reaction-diffusion model is suitable for the simulation of PLA degradation. In addition, simulation results by traditional reaction-diffusion model are also displayed in Figure 2 with dashed lines. It can be affirmed that traditional reaction-diffusion model is not appropriate to simulate PLA’s degradation with high mass loss.


12

Page 13 of 42

(b) 100

80

80

60

60

Mass loss (%)

Normalized molecular weight (%)

(a) 100

40

40 20

20 0

0

5

10

15

20

0

25

0

5

15

20

25

30

(d) 100

(c) 100

80

60

60

Mass loss (%)

80

40

40 20

20 0

10

Degradation time (week)


Normalized molecular weight (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

5

10

15

20

25

30

0

0

5

10

15

20

25

30



Figure 2. Normalized weight averaged molecular weights and mass losses as functions of time for the 2mm PLA plate ((a) and (b)) and the 0.3mm PLA film ((c) and (d)). The discrete symbols are experimental data 10, the solid lines and the dashed lines are fitted using the acceleratable reaction-diffusion model and traditional reaction-diffusion model

Table 1 Primary simulation parameters used in the studies


13


Acceleratable model

Page 14 of 42

Traditional model Autocataly

Parameters 0.3m

2mm

sis

Mixing

Noncatalysis

2mm

0.3mm

Dimension (no unit)

One

One

One

One

Two

Two

Two

M or M×N

2000

300

2000

300

100×100

100×100

100×100

Initial chains (103)

20000

1776

20000

1776

59200

59200

59200

Sample size (mm/mm2)

2

0.3

2

0.3

3×3

3×3

3×3

Time step Δt (hour)

1

1

1

1

1

1

1

Total time (week)

30

30

30

30

52

52

52

Calculation time (hour)

2

Cellular Automaton Simulation for Degradation of ... - ACS Publications

Recommend Documents