Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
On the initial propagation of dental caries Rene Fabregas and Jacob Rubinstein J. R. Soc. Interface 2014 11, 20140809, published 17 September 2014 References
This article cites 19 articles, 4 of which can be accessed free http://rsif.royalsocietypublishing.org/content/11/100/20140809.full#BIBL
Subject collections
Articles on similar topics can be found in the following collections: biomathematics (326 articles) biophysics (368 articles) computational biology (332 articles)
Email alerting service
Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here
To subscribe to J. R. Soc. Interface, go to: http://rsif.royalsocietypublishing.org/subscriptions
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
On the initial propagation of dental caries Rene Fabregas1,2 and Jacob Rubinstein1 1 2
Department of Mathematics, Technion-Israel Institute of Technology, Amado Building, Haifa 32000, Israel Departamento de Matema´tica Aplicada, Facultad de Matema´ticas, Universidad Complutense, Madrid 28040, Spain
rsif.royalsocietypublishing.org
Research Cite this article: Fabregas R, Rubinstein J. 2014 On the initial propagation of dental caries. J. R. Soc. Interface 11: 20140809. http://dx.doi.org/10.1098/rsif.2014.0809
Received: 22 July 2014 Accepted: 27 August 2014
Subject Areas: biomathematics, computational biology, biophysics Keywords: caries, enamel, mathematical model, prismless layer
Author for correspondence: Rene Fabregas e-mail: [email protected]
A multi-dimensional model for dental caries is applied to study the shape of caries lesions in a realistic tooth geometry and to examine the rate of progress of caries. An upgraded model, taking into account the outer prismless enamel layer, is derived and solved. The model demonstrates the importance of this layer in delaying the onset of caries. The conclusions are discussed in light of experimental results.
1. Introduction Caries formation is a fundamental process in dental health. It involves a complex array of chemical and transport mechanisms. As a first step towards the simulation of caries in a realistic geometry, and in order to develop tools for simulating treatments to arrest or prevent dental caries, we recently developed a model for this process [1]. Our derivation starts at the local transport and reaction processes. We averaged the equations over the local microscopic scale of the enamel prisms to obtain macroscopic equations that describe the concentrations of the main compounds and the dissolution of the enamel prisms. In particular, our model incorporates the anisotropic nature of the local geometry. We demonstrated the model in a simple two-dimensional rectangular geometry that is sometime used in experiments. This paper has two purposes. One is to apply the model to simulate caries progress in a more realistic tooth geometry, to identify the dependence of the lesion shape on the origin of caries on the tooth boundary and to examine carefully the question of caries propagation rate. A second goal is to incorporate into the model the prismless outer enamel layer and to examine its effect on caries progress. The dissolution of enamel follows the reaction of the hydroxyapatite molecules with hydrogen ions. The reaction is k1
2þ 2 Ca10 (PO4 )6 (OH)2(s) þ 8Hþ (aq) O 10Ca(aq) þ 6HPO4(aq) þ 2H2 O(l) : Hydroxyapatite (HAP)
k1
(1:1)
This reaction is reversible, but when the pH level drops to approximately 5 and lower, it proceeds mostly to the right [2]. As the enamel is reacting and is demineralized, material is lost and later cavities form in the enamel. Our model consists of equations for the evolution of the hydrogen and calcium ions, coupled with the evolution of the microstructure, i.e. the enamel prisms. The propagation of lesions was measured and analysed by a number of techniques. We draw particular attention to the work of Bjorndal and his co-worker [3,4]. The shape of the propagating lesion was summarized by Kidd & Fejerskov [5]. When the lesion starts at an approximal surface, it has a triangular shape in the enamel with a base at the enamel surface and a vertex at the enamel –dentine (ED) junction. Another ‘triangle’ is then observed at the dentine, with a base at the ED junction and a vertex deep in the dentine. On the other hand, when the lesion starts from a biofilm at an occlusal surface, i.e. a pit at the central top part of the tooth, it has the shape of a triangle with a vertex at the pit (enamel surface) and a base at the ED junction. Simulations of our model presented in this paper are in agreement with these observations. In §2, we shall present the model and the underlying assumptions behind it. In §3, we solve the equations in a three-dimensional rectangular domain. In §4, we solve the equations in a two-dimensional geometry that mimics a realistic
& 2014 The Author(s) Published by the Royal Society. All rights reserved.
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(b)
prisms
inter-prism region
rsif.royalsocietypublishing.org
E
Figure 1. (a) A sketch of realistic tooth geometry. The enamel and dentine regions are denoted by E and D, respectively. (b) A sketch of the prism array. tooth shape. In particular, we discuss there the important and controversial issue of what is the progression rate of caries. We also consider there the shape of the lesion, depending on where it starts along the tooth boundary. In §5, we upgrade the model to account for the resilient outer layer of the enamel and inspect its important effect on caries progress. Finally, we discuss the model and the simulations in §6 and propose future directions for it. We point out that there are several mathematical models for the initial stages of dental caries including references [6– 9] and so on. The models differ in the level of complexity (e.g. the number of compounds) and details of the reaction. These models are in general one-dimensional and thus anisotropy plays no role. Also, they are written down from the start on a macroscopic level.
2. The three-dimensional model A general sketch of the enamel and dentine parts of the tooth is depicted in figure 1a. The outer layer of the tooth consists of enamel in the form of a packed array of rods ( prisms) with radius of approximately 223 mm. The prisms vary considerably in length (2021000 mm, [10]), and the inter-prism distance is approximately 123 mm [10,11]. The prism long axis is in general orthogonal to the tooth surface. A sketch of the prisms is depicted in figure 1b. As we deal with a problem involving several length scales, we start with identifying them. Two length scales are geometrical. The first is l ¼ 1 mm, which corresponds to the prisms radius and separation between the prisms. A second geometrical length scale is the total thickness of the enamel layer, which is of the order of a few millimetres. In addition to these geometrical length scales, we identify a dynamical length scale, which is the scale where the physical processes taking place in the caries process are balanced, as follows. There are two main processes involved in caries progress. One is diffusion in the inter-prism space and the other one is reaction on the prism surface. Let D and R denote the local diffusion coefficient and the reaction rate at the prism boundary, respectively. As the dimensions of D and R are L2 =T and T , respectively, where L and T ffidenote length and time, pffiffiffiffiffiffiffiffiffi respectively, the term l ¼ D=R has the dimension of length. To understand the physical meaning of l, consider a simple one-dimensional model for diffusion and reaction in equilibrium, given by the equation 0 , x , 1,
u(0) ¼ u0 ,
u(1) ¼ 0,
where u is the concentration of some substance. The solution is u0e2x/l. Thus, l represents the size of the interval that is most affected by the variation of the boundary concentration. Using D ¼ 5 102 mm2 s21 and R ¼ 5 1022 l s21 [12,13], we obtain l 100 mm. The small ratio l/l implies a separation of scales. In what follows we shall use a nondimensional formulation, where length is scaled by l and time is scaled by 1/R. However, when we consider in §5 the effect of the prismless enamel outer layer, we use a much longer time scale as will be explained below. The macroscopic model consists of reaction and diffusion equations for two compounds: the calcium ion concentration c(x1, x2, x3, t) and the hydrogen ion concentration H(x1, x2, x3, t). In addition, we assume that the prisms essentially retain their cylindrical shape [14]; thus, the model includes an equation for the function r(x1, x2, x3, t), which is the local radius of a two-dimensional cross section of the prism. A detailed derivation of our model can be found in [1]. However, to make this paper self-contained, we explain the main features of the derivation and the meaning of the different terms in the equations. The first step is to characterize the local (microscopic) geometry. As described above, the prisms are roughly arranged in layers, where the long axis is perpendicular to the tooth surface. This implies that the global equations are inherently anisotropic, and, moreover, they depend on the shape of the tooth boundary. Some assumptions should be made on the prisms distribution. We assumed that the prisms form a locally periodic pattern in the planes orthogonal to their long axis. In the second step, we wrote down the local transport equations, including the reaction (1.1) that takes place at the surface of the prisms, and the diffusion of the different compounds in the interprism region. We next used the separation of scales, i.e. l l to perform a homogenization calculation [15,16]. The homogenization technique consists of introducing two different, essentially independent length scales, one for the micro-geometry and other for the macro-geometry, and then expanding the local equations in the small parameter l/l. The expansion splits the transport into a local one, which is integrated explicitly, and a global transport model, where the transport coefficients depend on the solution of the local equations. The problem of dental caries involves an unusual difficulty; the microscopic geometry is not fixed in time and space, as the prisms dissolve during demineralization. This means that the coefficients of the macroscopic equations depend on the evolution of the local microstructure. Thus, while the scales are separated, and so are the transport
J. R. Soc. Interface 11: 20140809
D
D@ 2x u(x) Ru(x) ¼ 0,
2
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
2 X
ij @ x )c þ j@ Ep jR1 (cs c) @ xi ( D j
i¼1
(2:1) and
@ t (jVjH) ¼ dH @ x3 (jVj@ x3 H) þ
2 X
! ij @ x )H @ xi ( D j
(2:2)
i¼1
j@ Ep jR1 (cs c): Here, @ xi and @ t denote partial differentiation, and dH is the ratio of the diffusion coefficients of the hydrogen and calcium ions. The reaction term is proportional to cs 2 c. The equilibrium concentration cs is modelled by cs ¼ gH:
(2:3)
The first term on the right-hand side in equations (2.1) and (2.2) models the diffusion along the ‘easy’ direction x3. The second term models the diffusion in the direction orthogonal can be to the prism long axis. The effective diffusion tensor D computed from the local geometry. Finally, the last term on the right-hand side of these equations models the averaged reaction taking place on the prism boundaries. The local geometry appears explicitly in equations (2.1)– (2.2) via the geometric term jVj, which is the local volume per unit length of the inter-prism domain, and j@ Epj, which is the local surface area per unit length of the prism. To fix ideas, and relying on the observations reported in [14], we assume that the prisms are essentially periodically arranged on a square lattice in the cross section orthogonal to the long axis, such that the prism centres are identified with the centres of the square lattice. All the local quantities are measured in units of l. For example, the lattice size scaled by l is 2. The local geometric terms are determined by r(x1, x2, x3, t), the local prism radius, which is also scaled by l, such that the maximal radius is 1. Note that when r ! 1, the lateral diffusion approaches zero. The geometric parameters in this setting are jVj ¼ 4 2 pr 2 and j@ Epj ¼ 2pr. The function r itself evolves according to the dissolution process
@ t r ¼ bR1 (gH c):
(2:4)
The model still contains three parameters: R1 is determined by the local reaction rate, g is determined by the precise model for the equilibrium concentration cs and b is related to the molar concentration of the hydroxyapatite crystals. However, we prefer to leave R1, g and b as free parameters that can later on be adjusted either from theoretical considerations or by fitting the model to observations.
3. Simulations in a three-dimensional rectangular geometry In this section, we consider a rectangular tooth sample. The domain is V ¼ {1 , x1 , x2 , 1, 0 , x3 , 3}: Before presenting specific simulations, we point out that the model above can be simplified if we approximate the factor dH by 1. Under this approximation, we can define the function w :¼ H c=g. The model then becomes
@ t (jVjw) ¼ @ x3 (jVj@ x3 w) þ
2 X
ij @ x )w (1 þ g)j@ Ep jR1 w, @ xi ( D j
i¼1
(3:1) together with the equation for the evolution of r
@ t r ¼ bgR1 w:
(3:2)
We point out that in reality the hydrogen diffusion coefficient is approximately five times larger than the calcium diffusion coefficient. However, numerical simulations of the full model give results that are qualitatively similar to the reduced model, and we believe that at the present level of caries modelling, the simplification gained by the approximation dH 1 is valuable. We assume the boundary condition w ¼ 0 at all the boundaries, except for the square x3 ¼ 0, (x1 , x2 ) [ [ 0:1, 0:1] [ 0:1, 0:1], which is the tooth outer boundary, where we apply the boundary condition w ¼ 1024. This boundary condition simulates a burst of low pH at the outer tooth surface, and the activity of the hydrogen ions is expected to start the caries process. We use the following parameters in the simulations: R1 ¼ 1, g ¼ 1 and b ¼ 3. The initial conditions are w ¼ 0 and r ¼ 0.9.
3
J. R. Soc. Interface 11: 20140809
@ t (jVjc) ¼ @ x3 (jVj@ x3 c) þ
As is readily seen from equations (2.1), the model is anisotropic, as the diffusion coefficient in the x3 direction (scaled to be 1 in our non-dimensional set-up) is different from the diffusion coefficients in the (x1, x2) directions, given by the effective This effective tensor D can be computed diffusion matrix D. separately from the transport equation on the microscopic scale [1]. As the microscopic geometry is fully encoded in the function r, this can be done in principle. However, practically, at any given point (x1, x2, x3, t) is a the problem of evaluating D formidable task. Fortunately, we found a way to greatly sim for many representative plify this task; we computed D values of r, and found out that, except for the limit of highly can be approximated by a linear function packed prisms, D of r: D(r) (4 3:64r)Id , where Id is the 2 2 identity matrix. We note that this simple formula does not hold when r is near 1, i.e. at highly packed prism configuration where the prisms nearly touch each other. At extreme values of r, a separate calculation should be done, solving in full the local [1]. Alternatively, asymptotic formula equations that govern D as those derived by Keller [17] can be used. We point out that modelling the progress of caries by following c, H and r enables direct correlation of the model with experiments. The importance of the hydrogen ions concentration H is evident. The concentration of the calcium ions c is useful because measuring the flux of the calcium ions leaving the tooth is a good indication of the amount of dissolved enamel. Another criterion for caries is the formation of observable lesions, and this can be correlated with the value of r.
rsif.royalsocietypublishing.org
processes, some coupling remains in the final model as we shall see below. The macroscopic model consists, therefore, of two equations for the evolution of the concentration of the calcium and hydrogen ions on the macroscopic (homogenized) scale. Although a comprehensive system of equations should be written down in general curvilinear coordinates that fit the global geometry, we find it more useful to write them down first for a rectangular domain, where the x3 axis is oriented along the prism long axis. Later, we shall formulate the model in a domain with spherical or circular symmetry that resembles a tooth shape. Thus, in a domain with rectangular symmetry, the model [1] is given by
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(c)
0.6
15.25
0.5
0.5
15.20
0.4
x3
0.3
0.4
0.2
0.1
0.1 500
1000
1500
0 15
15.15
E 15.10
0.3
0.2
0
~
4
15.05 15.00 20
t
25
30
35
40
500
4. Simulations in a circular geometry To simulate the model in a more realistic geometry, as sketched in figure 1a, we need in principle to construct a global system of curvilinear coordinates, as explained in §2, and write the model in these coordinates. However, the formulation can be greatly simplified if we consider that caries is, at least initially, a local process and moreover, caries typically starts either near the occlusal surface or near the approximal surface (dotted regions in figure 1a). In both cases, we can work locally in a spherical geometry (r, u, f ). Moreover, as the radial direction is orthogonal to the tooth surface in both cases, the ‘easy’ axis, where diffusion is faster, is along the radial variable r. Therefore, in such a geometry the model takes the form
@ t (jVjw) ¼
1 1 sin u @ u )w @ r (jVjr2 @ r w) þ 2 @ u (D r2 r sin u 1 @ f )w (1 þ g)j@ Ep jR1 w, þ 2 2 @ f (D r sin u
(4:1)
where again we use the simplified model with w ¼ H 2 c/g. As a further simplification, and also as another example of the model in a curvilinear geometry, we interpret the drawing in figure 1a as a two-dimensional sketch and solve the equation in a polar coordinate framework (r, u). The equation for w takes the form 1 1 @u )w (1 þ g)j@ Ep jR1 w: @ t (jVjw) ¼ @ r (rjVj@ r w) þ 2 @ u (D r r (4:2)
1500
pffi t and (c) the evolution
We solve equation (4.2), together with the equation for r and using the formula for jVj and j@ Epj in several set-ups. In all of them, the geometry is an angular sector G ¼ {(r, u): r1 , r , r2 , juj , b}:
(4:3)
We simulate caries progress when it starts at an occlusal surface. The boundary conditions are w ¼ 0 at all the boundaries, except at the sector Gs ¼ f(r,u)jr ¼ r1, juj , u0g, where w(r1, u, t) ¼ 1024. We use the parameters r1 ¼ 1, r2 ¼ 3, b ¼ p/4, u0 ¼ p/30, g ¼ 1, b ¼ 3 and R1 ¼ 1. Again, the goal is to simulate initial caries progress following a burst of low pH in a small portion of the boundary at an occlusal surface. Results: Just as in the previous example, we follow the motion x3(t) of the tip of a level set of r, as an indication of a visible lesion (figure 3a), as well as the total amount of ~ (figure 3d ). The tip’s evolution is also drawn as a enamel E pffiffi function of t to examine the importance of the diffusion process in the caries progress. We also draw in figure 3c the boundary of a given level set u(t). The shape of a caries lesion is often determined visually by inspecting a formation of a ‘white spot’, indicating a decrease in the total prism volume. As was pointed out in the Introduction, when a lesion starts near an occlusal surface, the lesion typically evolves in time into the enamel bulk by expanding its width, such that subsequent snapshots form a growing triangular shape. As the visible lesion reached the ED boundary, it expands in the horizontal direction. If we set arbitrarily the level set drawn in figure 3c to be the visual limit of the lesion, we obtain that the lesion boundary indeed propagates out by expanding from a small area at an occlusal surface. There is a debate in the literature about the actual rate of caries progression. Anderson et al. [18] cite a number of authors who argue that caries progresses as a linear function pffiffi of t. A typical result along this line is the paper of Poole et al. [19]. On the other hand, some authors, and in particular Anderson himself, provide experimental evidence that caries formation progresses linearly in t. The main argument in pffiffi favour of the t progression rate is that caries formation is dominated by the diffusion of the hydrogen ions, and as is well known, a particle undergoing Brownian motion propagates at this rate. However, this argument ignores the effect of the reaction term. Moreover, Stumpf & Porter [20] recently warned against a tendency in life sciences to associate observations with power laws. Indeed, the curves depicted in figures 2b and 3b are not convincing that the lesion propagates at a rate proportional
J. R. Soc. Interface 11: 20140809
Figure 2. (a) Motion of the tip of the fixed level set r ¼ 0.86, denoted by x3(t), (b) motion of tip of the same level set as a function of of the total enamel volume ~E(t).
Results: We present the progression of caries in two ways. First, we follow the tip of a certain level set of r. This mimics a visual inspection of the lesion. An alternative (but not equivalent) way is to follow the total amount of enamel. This is indicated by the amount of calcium ions leaving the tooth. The motion of the tip x3(t) of the level set r ¼ 0.86 is depicted in figure 2a. Recall that the basic time unit is 20 s, so the curve follows a period of about 10 h. The total ~ is amount of enamel, determined by the total volume E depicted in figure 2c. A number of authors argue that in light of the importance of the diffusion transport process, and as the distance covered by a Brownian particle is propffiffiffiffi portional to dt, where d is a diffusion coefficient, then pffiffi caries should propagate at rate proportional to t. Therefore, we depict in figure 2b the motion of the tip x3 as a function of pffiffi t. The curve seems not far from linear, but as we argue below this shape is misleading.
1000
t
t
rsif.royalsocietypublishing.org
x3
(b) 0.6
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(c)
0.6
0.5
0.5
0.2 0.1
0.4
0.4
r
q
0.3
0.3
0.2
0.2
0 –0.1
0.1
0.1
–0.2
0
500
1000
0 15
1500
t
20
25
30
35
40
500
1000
1500
t
–0.3
t
34.6
J. R. Soc. Interface 11: 20140809
(d )
5
0.3
rsif.royalsocietypublishing.org
r
(b)
0.6
34.5 34.4 ~
E
34.3 34.2 34.1 34.0 500
1000
1500
t
Figure 3. Progression of the tip of a level set r ¼ 0.86 of the enamel density function r: (a) motion of the tip of a level set, (b) motion of the same level set as a pffi function of t , (c) width of the level set 0.86 and (d ) the evolution of the total enamel volume ~E(t).
(a)
(b) 1.8
(c)
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1.7 1.6 1.5 r 1.4 1.3 1.2 1.1 20
30
40
50
60
70
4
5
t
6 t
7
8
20
30
40
50
60
t
Figure 4. Progression of the tip of a level set r1 ¼ 210: (a) motion of the tip of a level set, (b) motion of the same level set as a function of theoretical curve for the motion of the level set tip and (d ) the evolution of the total enamel volume ~E(t). pffiffi to t. Also, we point out that if we measure caries progression not by the propagation of a level set of the enamel density, but by the total amount of enamel removed from the tooth, then the progression looks linear in t (see figures 2c and 3d). To demonstrate the danger in such power laws, we simulated a case where the evolution of the tip of a level set of r can be solved analytically. We use the same model and geometry as in the last example, except that the boundary condition on the curve G1 (defined by r ¼ r1) is now w(r1, u, t) ¼ 1024 cos(2u). In light of the relatively slow enamel dissolution process, we solve the equation for w at a fixed geometry where r ¼ 0.98. At this extreme value of Instead, r, we cannot use anymore the linear formula for D. we evaluated the effective diffusion coefficient by solving the microscopic cell problem introduced in [1] and found ¼ 0:25. In addition, we used the parameters R1 ¼ 1/3, D g ¼ 3 and b ¼ 1/3. For this choice of parameters, the equations for w and r are approximately 1 1 @ t w ¼ @ r (r@ r w) þ 2 @ 2uu w 4w and @ t r ¼ w: r 4r
(4:4)
70
pffi t , (c) the
The function w reaches a steady state that can be found explicitly w ¼ 104
K1 (2r) cos(2u), K1 (2)
(4:5)
where K1 is the modified Bessel functions of the second kind. It is convenient to define r ¼ r0 þ 1024r1 where for the parameters we used @ tr1 ¼ 2104w. Substituting the expression given by w in equation (4.5) into the equation @ tr1 ¼ 2104w and solving for r1 at the tip, that is for u ¼ 0, we obtain
r1 (r, 0, t) ¼
K1 (2r) t: K1 (2)
(4:6)
It follows that the tip of the level set r1 ¼ 10 is given by the relation t(r) ¼ (K1(2)/K1(2r))10. Results: We use the formula for t(r) to draw the theoretical inverse curve r(t) in figure 4c. This curve can now be compared to the numerical solution r(t) given in figure 4a and pffiffi the curve r t depicted in figure 4b. It is clear from these pffiffi figures that the t is not a correct power law.
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(b)
1.0
6 34.5
rsif.royalsocietypublishing.org
0.2
0.8
34.0
0.1
r 0.6
~
E
q
0.4
400
600
t
–0.1
0.2 0
(c) 0.3
800
1000
33.0
–0.2 500
1000
t
1500
–0.3
(d)
33.5
500
1000
t
1500
0.2 0.1
q
400 –0.1
600
t
800
1000
–0.2 –0.3
Figure 5. Progression of the tip of a fixed level set 0.86 for the case of a lesion starting at an approximal surface of the tooth: (a) motion of the tip of a level set of r, (b) width of the level set 0.86, (c) the evolution of the total enamel volume ~E(t) and (d ) width of the level set 0.86 with fixed D ¼ 0:25. We examined above the progress of caries near an occlusal surface. We proceed to consider the shape of the caries lesion when it starts at an approximal surface. The domain is again the angular sector G; however, here the boundary conditions are w ¼ 0 at all the boundaries, except at the sector on G2 defined by Gs ¼ fr, u)jr ¼ r2, juj , u0g, where w(r2, u, t) ¼ 1024. We use the parameters r1 ¼ 1, r2 ¼ 3, b ¼ p/4, u0 ¼ p/30, bg ¼ 3 and R1 ¼ 1. Results: As in the previous examples, we follow the caries progress by drawing (figure 5a) the function x3(t), the tip of ~ which the level set r ¼ 0.86 and (figure 5c) the function E(t) is the total volume of enamel. In figure 5b, we draw the width of a given level set. An inverted triangular shape, as mentioned in §1 can be observed in this graph, although it is not as transparent as the triangular shape of a lesion starting at an occlusal surface (figure 3c). To understand why the triangular shapes for lesions starting at occlusal and approximal surfaces are different, it is useful to look more closely at the mechanism that underlies the progress of caries. When the lesion starts at an occlusal surface, the easy, that is radial, directions expand because of the negative curvature of the boundary. In addition, the lateral diffusion, although weak in light of the relatively smaller diffusion coefficient, also acts to expand the lesion’s width. Thus, both these mechanisms contribute to the triangular shape. On the other hand, when a lesion starts at an approximal surface, these mechanisms compete with each other. The easy, radial, directions now contract the lesion’s width as the boundary curvature is positive, while the lateral diffusion acts to expand the width. To illustrate this scenario, we repeated the last simulation with higher initial enamel density, and thus smaller lateral diffusion D. We follow the same level set and depict its width in figure 5d. Comparing figure 5b with 5d indeed shows that reducing the lateral diffusion enhances the inverted triangular shape.
5. The sound outer layer of enamel The outer layer of the enamel part of the tooth consists of a layer of hydroxyapatite molecules that are not organized in
prisms as in bulk enamel, but in a much denser and disordered pattern. In addition, this layer is known to be resistive to dissolution relative to bulk enamel. It is termed ‘relative sound enamel layer’ by Holly & Gray [6] and ‘prismless enamel layer’ by Ripa et al. [21] and by others. To incorporate this very important feature of tooth micro-geometry into our equations, and following the histological studies of Ripa et al. [21], we model this layer as a medium with a homogeneous and isotropic diffusion tensor, where the single diffusion coefficient is much smaller than the diffusion coefficients in bulk enamel. Also, as this outer layer is quite thin, and also more resistive to melting, we neglect the enamel dissolution in this layer. Therefore, we replace in this layer the reaction–diffusion equation (3.1) by the equation
@ t w ¼ r(Dr rw):
(5:1)
The thickness of the outer layer, that we denote here j, is reported to be in the range of 15–200 mm [6,21,22]. For the purpose of our simulation below, we take the layer to be of size 70 mm, which is approximately 0.7 in the non-dimensional length unit scaled by l. We could not find in the literature direct estimates for the diffusion coefficient Dr in the prismless layer. However, the experiments of Anderson et al. [18] give some indication for the magnitude of Dr. The authors report on a delay in the initiation of the enamel dissolution curve (fig. 4 in their paper). The delay is approximately 10–20 h. If we attribute this delay to the prismless outer layer, then using the estimate Dr j 2/T, and taking T ¼ 10 h and j ¼ 70 mm, we obtain Dr 1.5 1021 mm2 s21. To conform with the equations for bulk enamel that are written down in the non-dimensional units introduced at the beginning of §2, we denote the nondimensional diffusion coefficient by dr, which is the ratio between Dr and the main diffusion coefficient in bulk enamel that we took earlier to be approximately 5 102 mm2 s21; we thus obtain dr 1/3000. The model consists now of equation (5.1) in the prismless layer, and equation (3.1) or (4.1) in the bulk enamel. However, these equations involve very different diffusion coefficients. Therefore, it is useful to rescale the time variable
J. R. Soc. Interface 11: 20140809
0.3
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(b)
7 13.7
1.5
13.6 ~
E
2.0
rsif.royalsocietypublishing.org
x3 (t)
(a)
1.0
13.5
0.5
13.4 0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
t
Figure 6. The influence of the prismless layer on the evolution of a caries lesion for the boundary condition (5.5): (a) progression of the tip of the level set r ¼ 0.898 (solid curve) of models (5.2)–(5.4), and the tip of the same level set (dotted curve) for the models (3.1) and (3.2); (b) the evolution of the total enamel volume ~E(t) for the models (5.2)–(5.4) (solid line) and (3.1) and (3.2) (dotted line).
(a)
(b)
13.75 13.70
1.5
1.0
~
E
x3 (t)
13.65 13.60 13.55 13.50
0.5
13.45 0.1
0.2
0.3
0.4
0.5
t
0.1
0.2
0.3
0.4
0.5
t
Figure 7. Progression of dental caries with and without the influence of the outer layer for Example 6. (a) Progression of the tip of the level set r ¼ 0.898 (solid curve) of models (5.2)– (5.4), and the tip of the same level set (dotted curve) for the models (3.1) and (3.2); (b) the evolution of the total enamel volume ~E(t) for the models (5.2)– (5.4) (solid line) and (3.1) and (3.2) (dotted line).
in all the equations through t ¼ t0 /dr. Then, the left-hand side of equation (3.1) (or (4.1)) can be neglected. We write the augmented model that includes the prismless layer in a rectangular geometry, and to simplify the notation we omit the prime over the new time variable. The model thus consists of the following two equations:
@ t (jVjw) ¼
3 X
@ xi (jVj@ xi w),
0 , x3 , j
(5:2)
i¼1
and 0 ¼ @ x3 (jVj@ x3 w) þ
2 X
ij @ x w) (1 þ g)j@ Ep jR1 w, x3 . j: @ xi ( D j
i¼1
(5:3) Note that in the prismless layer, the geometry, i.e. V, is fixed; we introduced it into equation (5.2) in order to compare the diffusion coefficients of the two equations. The solution w as well as its flux are continuous across the interface x3 ¼ j. In addition to equations (5.2) and (5.3), we need to follow the evolution of r, although only in the bulk enamel x3 . j. In the new time scale, this equation reads
@tr ¼
bgR1 w: dr
(5:4)
Physically, the model means that the reaction –diffusion process in the bulk enamel is slaved by the boundary
behaviour of the prismless enamel, in the sense that it is always at equilibrium with the data of w at the boundary between the prismless strip and the bulk enamel. Also, while the basic time unit for the earlier model (3.1) was approximately 20 s, the time unit for the upgraded model (5.2) and (5.3) is approximately 16 h. Numerical solution of the system (5.2)– (5.4) is not trivial owing to the change of the nature of the equations at the boundary x3 ¼ j. We use here the method of domain decomposition with no overlapping (e.g. [23,24]). Two physiological set-ups were simulated using the model presented in this section. We use the same rectangular geometry as in the first example above. The prismless enamel domain is defined as x3 , j. We set as before equilibrium boundary condition w ¼ 0 at all boundaries except at the tooth surface x3 ¼ 0. There we set w(x1 , x3 ¼ 0, t) ¼ a(t)104 A(0:1 jx1 j),
(5:5)
where A is the Heaviside step function and a(t) ¼ 1 for 0 , t , 0.1, 0.5 , t , 0.7 and a(t) ¼ 0 for other values of t. This boundary condition mimics two localized bursts of low pH, say, following two meals separated by a dimensional time interval of approximately 6 h. Results: In figure 6, we draw the evolution of the tip of level set r ¼ 0.898. For comparison, we draw this level set
J. R. Soc. Interface 11: 20140809
t
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
A three-dimensional model for the initial evolution of dental caries was considered. We took into account the diffusion of two main compounds, hydrogen ions H and calcium ions c, and the reaction at the surface of the enamel prisms. The model consists of two transport equations (2.1) and (2.2) for c and H, together with equation (2.4) for the evolution of the enamel prisms on the microscopic level.
Funding statement. This work is supported in part by an F.I.R.S.T. Initial Training Grant of the ERC, and in part by a grant from the ISF.
References 1.
2.
3.
4.
5.
6.
Fabregas R, Rubinstein J. In press. A mathematical model for the progression of dental caries. Math. Med. Biol. (doi:10.1093/imammb/ dqt008) Shellis RP, Wilson RM. 2004 Apparent solubility distributions of hydroxyapatite and enamel apatite. J. Col. Interf. 278, 325 –332. (doi:10.1016/j.jcis. 2004.06.016) Bjorndal L, Thylstrup A. 1995 A structural analysis of approximal enamel caries lesions and subjacent dentin reactions. Eur. J. Oral Sci. 103, 25 –31. (doi:10.1111/j.1600-0722.1995. tb00006.x) Bjorndal L. 2008 The caries process and its effect on the pulp: the science is changing and so is our understanding. Pediatr. Dent. 7S, S2–S5. (doi:10. 1016/j.joen.2008.02.037) Kidd EAM, Fejerskov O. 2004 What constitutes dental caries? Histopathology of carious enamel and dentin related to the action of cariogenic biofilms. J. Dent. Res. 83, C35 –C38. (doi:10.1177/154405 910408301S07) Holly FJ, Gray JA. 1968 Mechanism for incipient carious lesion growth utilizing a physical model based on diffusion concept. Arch. Oral
Biol. 13, 319 –334. (doi:10.1016/0003-9969(68) 90130-1) 7. Patel MV, Fox JL, Higuchi WI. 1987 Physical model for non-steady-state dissolution of dental enamel. J. Dent. Res. 66, 1418 –1424. (doi:10.1177/ 00220345870660090201) 8. Patel MV, Fox JL, Higuchi WI. 1987 Effect of acid type on kinetics and dissolution of dental enamel demineralization. J. Dent. Res. 66, 1425–1430. (doi:10.1177/00220345870660090301) 9. Zimmerman SO. 1966 A mathematical theory of enamel solubility and the onset of dental caries: III. Development and computer simulation of a model of caries formation. Bull. Math. Biophys. 28, 443 –464. (doi:10.1007/BF02476825) 10. Radianski RJ, Renz H. 2004 A possible interdependency between the wavy path of enamel rods, distances of Retzius lines, and mitotic activity at the cervical loop in human teeth: a hypothesis. Med. Hypotheses 62, 945–949. (doi:10.1016/j. mehy.2003.12.038) 11. Jenng YR, Lin TT, Shieh DB. 2009 Nanotribological characterization of tooth enamel rod affected by surface treatment. J. Biomech. 42, 2249– 2254. (doi:10.1016/j.jbiomech.2009.06.057)
12. Bollet-Quivogne SEC, Anderson P, Dowker SEP, Eliott JC. 2005 Microradiographical study on the influence of diffusion in the external liquid on the rate of demineralization in hydroxyapatite aggregates. Eur. J. Oral Sci. 113, 53 –59. (doi:10.1111/j.16000722.2005.00188.x) 13. Van Dijk JWE, Borggreven JMPM, Driessens FCM. 1983 Diffusion in mammalian tooth enamel in relation to the caries process. Arch. Oral Biol. 28, 391–397. (doi:10.1016/0003-9969(83)90006-7) 14. Dumont ER. 1995 Mammalian enamel prism patterns and enamel deposition rates. Scanning Microsc. 9, 429–442. 15. Hornung U. 2012 Homogenization and porous media. Berlin, Germany: Springer. 16. Bensoussan A, Lions JL, Papanicolaou GC. 2011 Asymptotic analysis for periodic structures. New York, NY: Chelsea. 17. Keller JB. 1963 Conductivity of a medium containing a dense array of perfectly conducting spheres or cylinders or nonconducting cylinders. J. Appl. Phys. 34, 991–993. (doi:10.1063/1.1729580) 18. Anderson P, Levinkind M, Elliott JC. 1998 Scanning microradiographic studies of rates of in vitro demineralization in human and bovine dental
8
J. R. Soc. Interface 11: 20140809
6. Conclusion
The equations were first solved in a rectangular geometry. Then we presented two simulations in a circular (spherical) geometry that is adequate for most initial forms of caries. We emphasized two alternative ways of quantifying caries evolution: first by visual inspection, which is modelled by following a level set of the prism density distribution, and second by measuring the total volume of enamel. The simulations performed above in circular geometry model caries progress near an occlusal surface and near an approximal surface. In both cases, the model captures the well-known triangular shapes of the global caries lesion. However, we argued that this triangular shape should be more pronounced for lesions starting at an occlusal surface. The basic model was then upgraded to incorporate the very important effect of the outer prismless layer of enamel. The extended equations model this layer as a thin domain with isotropic diffusion matrix that is resistive to dissolution relating to bulk enamel. The presence of even a very thin such layer is shown to considerably reduce the progress of caries. A very interesting feature of the model is a qualitative explanation of fig. 4 of [18]. The model presented here still contains a few undetermined parameters, such as b, R1 and g. They can be estimated by matching the model to experiments. Then the model can be used to predict outcomes of new experiments, and in general, for in silico tests of different aspects of dental caries. There are a number of ways in which the model can be further extended. One is to improve the understanding of boundary conditions. Another very interesting direction is to add into the model the effect of fluoride.
rsif.royalsocietypublishing.org
(solid line) together with the evolution of the same level set when there is no prismless outer layer (dotted line). The powerful effect of the prismless layer and its role in slowing caries evolution is clearly seen in the figure. As a further demonstration of the important role played by the prismless outer layer, we apply a boundary condition that mimics a burst of localized low pH, say following a meal, followed by a period of elevated pH that gives rise to remineralization. The geometry is rectangular as before, with the prismless enamel domain defined as x3 , j. We set equilibrium boundary condition w ¼ 0 at all boundaries except at the tooth surface x3 ¼ 0. There we set w(x1, x3 ¼ 0, t) ¼ 1024A(0.1 2 jx1j) for 0 , t , 0.2, then w(x1, x3 ¼ 0, t) ¼ 21025 for 0.2 , t , 0.4 and finally, w(x1, x3 ¼ 0, t) ¼ 0 for 0.4 , t , 0.5 to complete the cycle of approximately 8 h. Results: In figure 7, we draw the evolution of the tip of level set r ¼ 0.898 (solid line) compared with the evolution of the same level set when there is no prismless outer layer (dotted line). The effect of the prismless layer is transparent. We also note that on the scale where the graph is drawn, there seems to be a delay in the initiation of the caries evolution owing to the outer layer.
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
21. Ripa LW, Gwinnnett AJ, Buonocore MG. 1966 The “prismless” outer layer of deciduous and permanent enamel. Arch. Oral Biol. 11, 41 –48. (doi:10.1016/ 0003-9969(66)90116-6) 22. Fava M, Watanabe I, Fava de Moraes F, Da Costa L. 1997 Prismless enamel in human non-erupted decidous moalr teeth: a scanning electron microscopy study. Rev. Odontol. Univ. Sao Paulo 11, 239–243. (doi:10.1590/S0103-06631997000400003)
23. Lions PL. 1990 On the Schwartz alternating method III: a variant for non-overlapping subdomains. In Third International symposium on domain decomposition methods for partial differential equations (eds TF Chan), pp. 202 –231. Phildelphia, PA: SIAM. 24. Quarteroni A, Valli A. 1999 Domain decomposition methods for partial differential equations. Oxford, UK: Oxford Science Publications.
9
rsif.royalsocietypublishing.org
enamel. Arch. Oral Biol. 43, 649– 656. (doi:10. 1016/S0003-9969(98)00052-1) 19. Poole DFG, Shellis RP, Tyler JE. 1981 Rates of formation in vitro of dental caries-like enamel lesion in man and some non-human primates. Arch. Oral Biol. 26, 413–417. (doi:10.1016/0003-9969(81)90039-X) 20. Stumpf MPH, Porter MA. 2012 Critical truths about power laws. Science 335, 665–666. (doi:10.1126/ science.1216142)
J. R. Soc. Interface 11: 20140809
On the initial propagation of dental caries Rene Fabregas and Jacob Rubinstein J. R. Soc. Interface 2014 11, 20140809, published 17 September 2014 References
This article cites 19 articles, 4 of which can be accessed free http://rsif.royalsocietypublishing.org/content/11/100/20140809.full#BIBL
Subject collections
Articles on similar topics can be found in the following collections: biomathematics (326 articles) biophysics (368 articles) computational biology (332 articles)
Email alerting service
Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here
To subscribe to J. R. Soc. Interface, go to: http://rsif.royalsocietypublishing.org/subscriptions
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
On the initial propagation of dental caries Rene Fabregas1,2 and Jacob Rubinstein1 1 2
Department of Mathematics, Technion-Israel Institute of Technology, Amado Building, Haifa 32000, Israel Departamento de Matema´tica Aplicada, Facultad de Matema´ticas, Universidad Complutense, Madrid 28040, Spain
rsif.royalsocietypublishing.org
Research Cite this article: Fabregas R, Rubinstein J. 2014 On the initial propagation of dental caries. J. R. Soc. Interface 11: 20140809. http://dx.doi.org/10.1098/rsif.2014.0809
Received: 22 July 2014 Accepted: 27 August 2014
Subject Areas: biomathematics, computational biology, biophysics Keywords: caries, enamel, mathematical model, prismless layer
Author for correspondence: Rene Fabregas e-mail: [email protected]
A multi-dimensional model for dental caries is applied to study the shape of caries lesions in a realistic tooth geometry and to examine the rate of progress of caries. An upgraded model, taking into account the outer prismless enamel layer, is derived and solved. The model demonstrates the importance of this layer in delaying the onset of caries. The conclusions are discussed in light of experimental results.
1. Introduction Caries formation is a fundamental process in dental health. It involves a complex array of chemical and transport mechanisms. As a first step towards the simulation of caries in a realistic geometry, and in order to develop tools for simulating treatments to arrest or prevent dental caries, we recently developed a model for this process [1]. Our derivation starts at the local transport and reaction processes. We averaged the equations over the local microscopic scale of the enamel prisms to obtain macroscopic equations that describe the concentrations of the main compounds and the dissolution of the enamel prisms. In particular, our model incorporates the anisotropic nature of the local geometry. We demonstrated the model in a simple two-dimensional rectangular geometry that is sometime used in experiments. This paper has two purposes. One is to apply the model to simulate caries progress in a more realistic tooth geometry, to identify the dependence of the lesion shape on the origin of caries on the tooth boundary and to examine carefully the question of caries propagation rate. A second goal is to incorporate into the model the prismless outer enamel layer and to examine its effect on caries progress. The dissolution of enamel follows the reaction of the hydroxyapatite molecules with hydrogen ions. The reaction is k1
2þ 2 Ca10 (PO4 )6 (OH)2(s) þ 8Hþ (aq) O 10Ca(aq) þ 6HPO4(aq) þ 2H2 O(l) : Hydroxyapatite (HAP)
k1
(1:1)
This reaction is reversible, but when the pH level drops to approximately 5 and lower, it proceeds mostly to the right [2]. As the enamel is reacting and is demineralized, material is lost and later cavities form in the enamel. Our model consists of equations for the evolution of the hydrogen and calcium ions, coupled with the evolution of the microstructure, i.e. the enamel prisms. The propagation of lesions was measured and analysed by a number of techniques. We draw particular attention to the work of Bjorndal and his co-worker [3,4]. The shape of the propagating lesion was summarized by Kidd & Fejerskov [5]. When the lesion starts at an approximal surface, it has a triangular shape in the enamel with a base at the enamel surface and a vertex at the enamel –dentine (ED) junction. Another ‘triangle’ is then observed at the dentine, with a base at the ED junction and a vertex deep in the dentine. On the other hand, when the lesion starts from a biofilm at an occlusal surface, i.e. a pit at the central top part of the tooth, it has the shape of a triangle with a vertex at the pit (enamel surface) and a base at the ED junction. Simulations of our model presented in this paper are in agreement with these observations. In §2, we shall present the model and the underlying assumptions behind it. In §3, we solve the equations in a three-dimensional rectangular domain. In §4, we solve the equations in a two-dimensional geometry that mimics a realistic
& 2014 The Author(s) Published by the Royal Society. All rights reserved.
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(b)
prisms
inter-prism region
rsif.royalsocietypublishing.org
E
Figure 1. (a) A sketch of realistic tooth geometry. The enamel and dentine regions are denoted by E and D, respectively. (b) A sketch of the prism array. tooth shape. In particular, we discuss there the important and controversial issue of what is the progression rate of caries. We also consider there the shape of the lesion, depending on where it starts along the tooth boundary. In §5, we upgrade the model to account for the resilient outer layer of the enamel and inspect its important effect on caries progress. Finally, we discuss the model and the simulations in §6 and propose future directions for it. We point out that there are several mathematical models for the initial stages of dental caries including references [6– 9] and so on. The models differ in the level of complexity (e.g. the number of compounds) and details of the reaction. These models are in general one-dimensional and thus anisotropy plays no role. Also, they are written down from the start on a macroscopic level.
2. The three-dimensional model A general sketch of the enamel and dentine parts of the tooth is depicted in figure 1a. The outer layer of the tooth consists of enamel in the form of a packed array of rods ( prisms) with radius of approximately 223 mm. The prisms vary considerably in length (2021000 mm, [10]), and the inter-prism distance is approximately 123 mm [10,11]. The prism long axis is in general orthogonal to the tooth surface. A sketch of the prisms is depicted in figure 1b. As we deal with a problem involving several length scales, we start with identifying them. Two length scales are geometrical. The first is l ¼ 1 mm, which corresponds to the prisms radius and separation between the prisms. A second geometrical length scale is the total thickness of the enamel layer, which is of the order of a few millimetres. In addition to these geometrical length scales, we identify a dynamical length scale, which is the scale where the physical processes taking place in the caries process are balanced, as follows. There are two main processes involved in caries progress. One is diffusion in the inter-prism space and the other one is reaction on the prism surface. Let D and R denote the local diffusion coefficient and the reaction rate at the prism boundary, respectively. As the dimensions of D and R are L2 =T and T , respectively, where L and T ffidenote length and time, pffiffiffiffiffiffiffiffiffi respectively, the term l ¼ D=R has the dimension of length. To understand the physical meaning of l, consider a simple one-dimensional model for diffusion and reaction in equilibrium, given by the equation 0 , x , 1,
u(0) ¼ u0 ,
u(1) ¼ 0,
where u is the concentration of some substance. The solution is u0e2x/l. Thus, l represents the size of the interval that is most affected by the variation of the boundary concentration. Using D ¼ 5 102 mm2 s21 and R ¼ 5 1022 l s21 [12,13], we obtain l 100 mm. The small ratio l/l implies a separation of scales. In what follows we shall use a nondimensional formulation, where length is scaled by l and time is scaled by 1/R. However, when we consider in §5 the effect of the prismless enamel outer layer, we use a much longer time scale as will be explained below. The macroscopic model consists of reaction and diffusion equations for two compounds: the calcium ion concentration c(x1, x2, x3, t) and the hydrogen ion concentration H(x1, x2, x3, t). In addition, we assume that the prisms essentially retain their cylindrical shape [14]; thus, the model includes an equation for the function r(x1, x2, x3, t), which is the local radius of a two-dimensional cross section of the prism. A detailed derivation of our model can be found in [1]. However, to make this paper self-contained, we explain the main features of the derivation and the meaning of the different terms in the equations. The first step is to characterize the local (microscopic) geometry. As described above, the prisms are roughly arranged in layers, where the long axis is perpendicular to the tooth surface. This implies that the global equations are inherently anisotropic, and, moreover, they depend on the shape of the tooth boundary. Some assumptions should be made on the prisms distribution. We assumed that the prisms form a locally periodic pattern in the planes orthogonal to their long axis. In the second step, we wrote down the local transport equations, including the reaction (1.1) that takes place at the surface of the prisms, and the diffusion of the different compounds in the interprism region. We next used the separation of scales, i.e. l l to perform a homogenization calculation [15,16]. The homogenization technique consists of introducing two different, essentially independent length scales, one for the micro-geometry and other for the macro-geometry, and then expanding the local equations in the small parameter l/l. The expansion splits the transport into a local one, which is integrated explicitly, and a global transport model, where the transport coefficients depend on the solution of the local equations. The problem of dental caries involves an unusual difficulty; the microscopic geometry is not fixed in time and space, as the prisms dissolve during demineralization. This means that the coefficients of the macroscopic equations depend on the evolution of the local microstructure. Thus, while the scales are separated, and so are the transport
J. R. Soc. Interface 11: 20140809
D
D@ 2x u(x) Ru(x) ¼ 0,
2
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
2 X
ij @ x )c þ j@ Ep jR1 (cs c) @ xi ( D j
i¼1
(2:1) and
@ t (jVjH) ¼ dH @ x3 (jVj@ x3 H) þ
2 X
! ij @ x )H @ xi ( D j
(2:2)
i¼1
j@ Ep jR1 (cs c): Here, @ xi and @ t denote partial differentiation, and dH is the ratio of the diffusion coefficients of the hydrogen and calcium ions. The reaction term is proportional to cs 2 c. The equilibrium concentration cs is modelled by cs ¼ gH:
(2:3)
The first term on the right-hand side in equations (2.1) and (2.2) models the diffusion along the ‘easy’ direction x3. The second term models the diffusion in the direction orthogonal can be to the prism long axis. The effective diffusion tensor D computed from the local geometry. Finally, the last term on the right-hand side of these equations models the averaged reaction taking place on the prism boundaries. The local geometry appears explicitly in equations (2.1)– (2.2) via the geometric term jVj, which is the local volume per unit length of the inter-prism domain, and j@ Epj, which is the local surface area per unit length of the prism. To fix ideas, and relying on the observations reported in [14], we assume that the prisms are essentially periodically arranged on a square lattice in the cross section orthogonal to the long axis, such that the prism centres are identified with the centres of the square lattice. All the local quantities are measured in units of l. For example, the lattice size scaled by l is 2. The local geometric terms are determined by r(x1, x2, x3, t), the local prism radius, which is also scaled by l, such that the maximal radius is 1. Note that when r ! 1, the lateral diffusion approaches zero. The geometric parameters in this setting are jVj ¼ 4 2 pr 2 and j@ Epj ¼ 2pr. The function r itself evolves according to the dissolution process
@ t r ¼ bR1 (gH c):
(2:4)
The model still contains three parameters: R1 is determined by the local reaction rate, g is determined by the precise model for the equilibrium concentration cs and b is related to the molar concentration of the hydroxyapatite crystals. However, we prefer to leave R1, g and b as free parameters that can later on be adjusted either from theoretical considerations or by fitting the model to observations.
3. Simulations in a three-dimensional rectangular geometry In this section, we consider a rectangular tooth sample. The domain is V ¼ {1 , x1 , x2 , 1, 0 , x3 , 3}: Before presenting specific simulations, we point out that the model above can be simplified if we approximate the factor dH by 1. Under this approximation, we can define the function w :¼ H c=g. The model then becomes
@ t (jVjw) ¼ @ x3 (jVj@ x3 w) þ
2 X
ij @ x )w (1 þ g)j@ Ep jR1 w, @ xi ( D j
i¼1
(3:1) together with the equation for the evolution of r
@ t r ¼ bgR1 w:
(3:2)
We point out that in reality the hydrogen diffusion coefficient is approximately five times larger than the calcium diffusion coefficient. However, numerical simulations of the full model give results that are qualitatively similar to the reduced model, and we believe that at the present level of caries modelling, the simplification gained by the approximation dH 1 is valuable. We assume the boundary condition w ¼ 0 at all the boundaries, except for the square x3 ¼ 0, (x1 , x2 ) [ [ 0:1, 0:1] [ 0:1, 0:1], which is the tooth outer boundary, where we apply the boundary condition w ¼ 1024. This boundary condition simulates a burst of low pH at the outer tooth surface, and the activity of the hydrogen ions is expected to start the caries process. We use the following parameters in the simulations: R1 ¼ 1, g ¼ 1 and b ¼ 3. The initial conditions are w ¼ 0 and r ¼ 0.9.
3
J. R. Soc. Interface 11: 20140809
@ t (jVjc) ¼ @ x3 (jVj@ x3 c) þ
As is readily seen from equations (2.1), the model is anisotropic, as the diffusion coefficient in the x3 direction (scaled to be 1 in our non-dimensional set-up) is different from the diffusion coefficients in the (x1, x2) directions, given by the effective This effective tensor D can be computed diffusion matrix D. separately from the transport equation on the microscopic scale [1]. As the microscopic geometry is fully encoded in the function r, this can be done in principle. However, practically, at any given point (x1, x2, x3, t) is a the problem of evaluating D formidable task. Fortunately, we found a way to greatly sim for many representative plify this task; we computed D values of r, and found out that, except for the limit of highly can be approximated by a linear function packed prisms, D of r: D(r) (4 3:64r)Id , where Id is the 2 2 identity matrix. We note that this simple formula does not hold when r is near 1, i.e. at highly packed prism configuration where the prisms nearly touch each other. At extreme values of r, a separate calculation should be done, solving in full the local [1]. Alternatively, asymptotic formula equations that govern D as those derived by Keller [17] can be used. We point out that modelling the progress of caries by following c, H and r enables direct correlation of the model with experiments. The importance of the hydrogen ions concentration H is evident. The concentration of the calcium ions c is useful because measuring the flux of the calcium ions leaving the tooth is a good indication of the amount of dissolved enamel. Another criterion for caries is the formation of observable lesions, and this can be correlated with the value of r.
rsif.royalsocietypublishing.org
processes, some coupling remains in the final model as we shall see below. The macroscopic model consists, therefore, of two equations for the evolution of the concentration of the calcium and hydrogen ions on the macroscopic (homogenized) scale. Although a comprehensive system of equations should be written down in general curvilinear coordinates that fit the global geometry, we find it more useful to write them down first for a rectangular domain, where the x3 axis is oriented along the prism long axis. Later, we shall formulate the model in a domain with spherical or circular symmetry that resembles a tooth shape. Thus, in a domain with rectangular symmetry, the model [1] is given by
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(c)
0.6
15.25
0.5
0.5
15.20
0.4
x3
0.3
0.4
0.2
0.1
0.1 500
1000
1500
0 15
15.15
E 15.10
0.3
0.2
0
~
4
15.05 15.00 20
t
25
30
35
40
500
4. Simulations in a circular geometry To simulate the model in a more realistic geometry, as sketched in figure 1a, we need in principle to construct a global system of curvilinear coordinates, as explained in §2, and write the model in these coordinates. However, the formulation can be greatly simplified if we consider that caries is, at least initially, a local process and moreover, caries typically starts either near the occlusal surface or near the approximal surface (dotted regions in figure 1a). In both cases, we can work locally in a spherical geometry (r, u, f ). Moreover, as the radial direction is orthogonal to the tooth surface in both cases, the ‘easy’ axis, where diffusion is faster, is along the radial variable r. Therefore, in such a geometry the model takes the form
@ t (jVjw) ¼
1 1 sin u @ u )w @ r (jVjr2 @ r w) þ 2 @ u (D r2 r sin u 1 @ f )w (1 þ g)j@ Ep jR1 w, þ 2 2 @ f (D r sin u
(4:1)
where again we use the simplified model with w ¼ H 2 c/g. As a further simplification, and also as another example of the model in a curvilinear geometry, we interpret the drawing in figure 1a as a two-dimensional sketch and solve the equation in a polar coordinate framework (r, u). The equation for w takes the form 1 1 @u )w (1 þ g)j@ Ep jR1 w: @ t (jVjw) ¼ @ r (rjVj@ r w) þ 2 @ u (D r r (4:2)
1500
pffi t and (c) the evolution
We solve equation (4.2), together with the equation for r and using the formula for jVj and j@ Epj in several set-ups. In all of them, the geometry is an angular sector G ¼ {(r, u): r1 , r , r2 , juj , b}:
(4:3)
We simulate caries progress when it starts at an occlusal surface. The boundary conditions are w ¼ 0 at all the boundaries, except at the sector Gs ¼ f(r,u)jr ¼ r1, juj , u0g, where w(r1, u, t) ¼ 1024. We use the parameters r1 ¼ 1, r2 ¼ 3, b ¼ p/4, u0 ¼ p/30, g ¼ 1, b ¼ 3 and R1 ¼ 1. Again, the goal is to simulate initial caries progress following a burst of low pH in a small portion of the boundary at an occlusal surface. Results: Just as in the previous example, we follow the motion x3(t) of the tip of a level set of r, as an indication of a visible lesion (figure 3a), as well as the total amount of ~ (figure 3d ). The tip’s evolution is also drawn as a enamel E pffiffi function of t to examine the importance of the diffusion process in the caries progress. We also draw in figure 3c the boundary of a given level set u(t). The shape of a caries lesion is often determined visually by inspecting a formation of a ‘white spot’, indicating a decrease in the total prism volume. As was pointed out in the Introduction, when a lesion starts near an occlusal surface, the lesion typically evolves in time into the enamel bulk by expanding its width, such that subsequent snapshots form a growing triangular shape. As the visible lesion reached the ED boundary, it expands in the horizontal direction. If we set arbitrarily the level set drawn in figure 3c to be the visual limit of the lesion, we obtain that the lesion boundary indeed propagates out by expanding from a small area at an occlusal surface. There is a debate in the literature about the actual rate of caries progression. Anderson et al. [18] cite a number of authors who argue that caries progresses as a linear function pffiffi of t. A typical result along this line is the paper of Poole et al. [19]. On the other hand, some authors, and in particular Anderson himself, provide experimental evidence that caries formation progresses linearly in t. The main argument in pffiffi favour of the t progression rate is that caries formation is dominated by the diffusion of the hydrogen ions, and as is well known, a particle undergoing Brownian motion propagates at this rate. However, this argument ignores the effect of the reaction term. Moreover, Stumpf & Porter [20] recently warned against a tendency in life sciences to associate observations with power laws. Indeed, the curves depicted in figures 2b and 3b are not convincing that the lesion propagates at a rate proportional
J. R. Soc. Interface 11: 20140809
Figure 2. (a) Motion of the tip of the fixed level set r ¼ 0.86, denoted by x3(t), (b) motion of tip of the same level set as a function of of the total enamel volume ~E(t).
Results: We present the progression of caries in two ways. First, we follow the tip of a certain level set of r. This mimics a visual inspection of the lesion. An alternative (but not equivalent) way is to follow the total amount of enamel. This is indicated by the amount of calcium ions leaving the tooth. The motion of the tip x3(t) of the level set r ¼ 0.86 is depicted in figure 2a. Recall that the basic time unit is 20 s, so the curve follows a period of about 10 h. The total ~ is amount of enamel, determined by the total volume E depicted in figure 2c. A number of authors argue that in light of the importance of the diffusion transport process, and as the distance covered by a Brownian particle is propffiffiffiffi portional to dt, where d is a diffusion coefficient, then pffiffi caries should propagate at rate proportional to t. Therefore, we depict in figure 2b the motion of the tip x3 as a function of pffiffi t. The curve seems not far from linear, but as we argue below this shape is misleading.
1000
t
t
rsif.royalsocietypublishing.org
x3
(b) 0.6
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(c)
0.6
0.5
0.5
0.2 0.1
0.4
0.4
r
q
0.3
0.3
0.2
0.2
0 –0.1
0.1
0.1
–0.2
0
500
1000
0 15
1500
t
20
25
30
35
40
500
1000
1500
t
–0.3
t
34.6
J. R. Soc. Interface 11: 20140809
(d )
5
0.3
rsif.royalsocietypublishing.org
r
(b)
0.6
34.5 34.4 ~
E
34.3 34.2 34.1 34.0 500
1000
1500
t
Figure 3. Progression of the tip of a level set r ¼ 0.86 of the enamel density function r: (a) motion of the tip of a level set, (b) motion of the same level set as a pffi function of t , (c) width of the level set 0.86 and (d ) the evolution of the total enamel volume ~E(t).
(a)
(b) 1.8
(c)
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1.7 1.6 1.5 r 1.4 1.3 1.2 1.1 20
30
40
50
60
70
4
5
t
6 t
7
8
20
30
40
50
60
t
Figure 4. Progression of the tip of a level set r1 ¼ 210: (a) motion of the tip of a level set, (b) motion of the same level set as a function of theoretical curve for the motion of the level set tip and (d ) the evolution of the total enamel volume ~E(t). pffiffi to t. Also, we point out that if we measure caries progression not by the propagation of a level set of the enamel density, but by the total amount of enamel removed from the tooth, then the progression looks linear in t (see figures 2c and 3d). To demonstrate the danger in such power laws, we simulated a case where the evolution of the tip of a level set of r can be solved analytically. We use the same model and geometry as in the last example, except that the boundary condition on the curve G1 (defined by r ¼ r1) is now w(r1, u, t) ¼ 1024 cos(2u). In light of the relatively slow enamel dissolution process, we solve the equation for w at a fixed geometry where r ¼ 0.98. At this extreme value of Instead, r, we cannot use anymore the linear formula for D. we evaluated the effective diffusion coefficient by solving the microscopic cell problem introduced in [1] and found ¼ 0:25. In addition, we used the parameters R1 ¼ 1/3, D g ¼ 3 and b ¼ 1/3. For this choice of parameters, the equations for w and r are approximately 1 1 @ t w ¼ @ r (r@ r w) þ 2 @ 2uu w 4w and @ t r ¼ w: r 4r
(4:4)
70
pffi t , (c) the
The function w reaches a steady state that can be found explicitly w ¼ 104
K1 (2r) cos(2u), K1 (2)
(4:5)
where K1 is the modified Bessel functions of the second kind. It is convenient to define r ¼ r0 þ 1024r1 where for the parameters we used @ tr1 ¼ 2104w. Substituting the expression given by w in equation (4.5) into the equation @ tr1 ¼ 2104w and solving for r1 at the tip, that is for u ¼ 0, we obtain
r1 (r, 0, t) ¼
K1 (2r) t: K1 (2)
(4:6)
It follows that the tip of the level set r1 ¼ 10 is given by the relation t(r) ¼ (K1(2)/K1(2r))10. Results: We use the formula for t(r) to draw the theoretical inverse curve r(t) in figure 4c. This curve can now be compared to the numerical solution r(t) given in figure 4a and pffiffi the curve r t depicted in figure 4b. It is clear from these pffiffi figures that the t is not a correct power law.
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(a)
(b)
1.0
6 34.5
rsif.royalsocietypublishing.org
0.2
0.8
34.0
0.1
r 0.6
~
E
q
0.4
400
600
t
–0.1
0.2 0
(c) 0.3
800
1000
33.0
–0.2 500
1000
t
1500
–0.3
(d)
33.5
500
1000
t
1500
0.2 0.1
q
400 –0.1
600
t
800
1000
–0.2 –0.3
Figure 5. Progression of the tip of a fixed level set 0.86 for the case of a lesion starting at an approximal surface of the tooth: (a) motion of the tip of a level set of r, (b) width of the level set 0.86, (c) the evolution of the total enamel volume ~E(t) and (d ) width of the level set 0.86 with fixed D ¼ 0:25. We examined above the progress of caries near an occlusal surface. We proceed to consider the shape of the caries lesion when it starts at an approximal surface. The domain is again the angular sector G; however, here the boundary conditions are w ¼ 0 at all the boundaries, except at the sector on G2 defined by Gs ¼ fr, u)jr ¼ r2, juj , u0g, where w(r2, u, t) ¼ 1024. We use the parameters r1 ¼ 1, r2 ¼ 3, b ¼ p/4, u0 ¼ p/30, bg ¼ 3 and R1 ¼ 1. Results: As in the previous examples, we follow the caries progress by drawing (figure 5a) the function x3(t), the tip of ~ which the level set r ¼ 0.86 and (figure 5c) the function E(t) is the total volume of enamel. In figure 5b, we draw the width of a given level set. An inverted triangular shape, as mentioned in §1 can be observed in this graph, although it is not as transparent as the triangular shape of a lesion starting at an occlusal surface (figure 3c). To understand why the triangular shapes for lesions starting at occlusal and approximal surfaces are different, it is useful to look more closely at the mechanism that underlies the progress of caries. When the lesion starts at an occlusal surface, the easy, that is radial, directions expand because of the negative curvature of the boundary. In addition, the lateral diffusion, although weak in light of the relatively smaller diffusion coefficient, also acts to expand the lesion’s width. Thus, both these mechanisms contribute to the triangular shape. On the other hand, when a lesion starts at an approximal surface, these mechanisms compete with each other. The easy, radial, directions now contract the lesion’s width as the boundary curvature is positive, while the lateral diffusion acts to expand the width. To illustrate this scenario, we repeated the last simulation with higher initial enamel density, and thus smaller lateral diffusion D. We follow the same level set and depict its width in figure 5d. Comparing figure 5b with 5d indeed shows that reducing the lateral diffusion enhances the inverted triangular shape.
5. The sound outer layer of enamel The outer layer of the enamel part of the tooth consists of a layer of hydroxyapatite molecules that are not organized in
prisms as in bulk enamel, but in a much denser and disordered pattern. In addition, this layer is known to be resistive to dissolution relative to bulk enamel. It is termed ‘relative sound enamel layer’ by Holly & Gray [6] and ‘prismless enamel layer’ by Ripa et al. [21] and by others. To incorporate this very important feature of tooth micro-geometry into our equations, and following the histological studies of Ripa et al. [21], we model this layer as a medium with a homogeneous and isotropic diffusion tensor, where the single diffusion coefficient is much smaller than the diffusion coefficients in bulk enamel. Also, as this outer layer is quite thin, and also more resistive to melting, we neglect the enamel dissolution in this layer. Therefore, we replace in this layer the reaction–diffusion equation (3.1) by the equation
@ t w ¼ r(Dr rw):
(5:1)
The thickness of the outer layer, that we denote here j, is reported to be in the range of 15–200 mm [6,21,22]. For the purpose of our simulation below, we take the layer to be of size 70 mm, which is approximately 0.7 in the non-dimensional length unit scaled by l. We could not find in the literature direct estimates for the diffusion coefficient Dr in the prismless layer. However, the experiments of Anderson et al. [18] give some indication for the magnitude of Dr. The authors report on a delay in the initiation of the enamel dissolution curve (fig. 4 in their paper). The delay is approximately 10–20 h. If we attribute this delay to the prismless outer layer, then using the estimate Dr j 2/T, and taking T ¼ 10 h and j ¼ 70 mm, we obtain Dr 1.5 1021 mm2 s21. To conform with the equations for bulk enamel that are written down in the non-dimensional units introduced at the beginning of §2, we denote the nondimensional diffusion coefficient by dr, which is the ratio between Dr and the main diffusion coefficient in bulk enamel that we took earlier to be approximately 5 102 mm2 s21; we thus obtain dr 1/3000. The model consists now of equation (5.1) in the prismless layer, and equation (3.1) or (4.1) in the bulk enamel. However, these equations involve very different diffusion coefficients. Therefore, it is useful to rescale the time variable
J. R. Soc. Interface 11: 20140809
0.3
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
(b)
7 13.7
1.5
13.6 ~
E
2.0
rsif.royalsocietypublishing.org
x3 (t)
(a)
1.0
13.5
0.5
13.4 0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
t
Figure 6. The influence of the prismless layer on the evolution of a caries lesion for the boundary condition (5.5): (a) progression of the tip of the level set r ¼ 0.898 (solid curve) of models (5.2)–(5.4), and the tip of the same level set (dotted curve) for the models (3.1) and (3.2); (b) the evolution of the total enamel volume ~E(t) for the models (5.2)–(5.4) (solid line) and (3.1) and (3.2) (dotted line).
(a)
(b)
13.75 13.70
1.5
1.0
~
E
x3 (t)
13.65 13.60 13.55 13.50
0.5
13.45 0.1
0.2
0.3
0.4
0.5
t
0.1
0.2
0.3
0.4
0.5
t
Figure 7. Progression of dental caries with and without the influence of the outer layer for Example 6. (a) Progression of the tip of the level set r ¼ 0.898 (solid curve) of models (5.2)– (5.4), and the tip of the same level set (dotted curve) for the models (3.1) and (3.2); (b) the evolution of the total enamel volume ~E(t) for the models (5.2)– (5.4) (solid line) and (3.1) and (3.2) (dotted line).
in all the equations through t ¼ t0 /dr. Then, the left-hand side of equation (3.1) (or (4.1)) can be neglected. We write the augmented model that includes the prismless layer in a rectangular geometry, and to simplify the notation we omit the prime over the new time variable. The model thus consists of the following two equations:
@ t (jVjw) ¼
3 X
@ xi (jVj@ xi w),
0 , x3 , j
(5:2)
i¼1
and 0 ¼ @ x3 (jVj@ x3 w) þ
2 X
ij @ x w) (1 þ g)j@ Ep jR1 w, x3 . j: @ xi ( D j
i¼1
(5:3) Note that in the prismless layer, the geometry, i.e. V, is fixed; we introduced it into equation (5.2) in order to compare the diffusion coefficients of the two equations. The solution w as well as its flux are continuous across the interface x3 ¼ j. In addition to equations (5.2) and (5.3), we need to follow the evolution of r, although only in the bulk enamel x3 . j. In the new time scale, this equation reads
@tr ¼
bgR1 w: dr
(5:4)
Physically, the model means that the reaction –diffusion process in the bulk enamel is slaved by the boundary
behaviour of the prismless enamel, in the sense that it is always at equilibrium with the data of w at the boundary between the prismless strip and the bulk enamel. Also, while the basic time unit for the earlier model (3.1) was approximately 20 s, the time unit for the upgraded model (5.2) and (5.3) is approximately 16 h. Numerical solution of the system (5.2)– (5.4) is not trivial owing to the change of the nature of the equations at the boundary x3 ¼ j. We use here the method of domain decomposition with no overlapping (e.g. [23,24]). Two physiological set-ups were simulated using the model presented in this section. We use the same rectangular geometry as in the first example above. The prismless enamel domain is defined as x3 , j. We set as before equilibrium boundary condition w ¼ 0 at all boundaries except at the tooth surface x3 ¼ 0. There we set w(x1 , x3 ¼ 0, t) ¼ a(t)104 A(0:1 jx1 j),
(5:5)
where A is the Heaviside step function and a(t) ¼ 1 for 0 , t , 0.1, 0.5 , t , 0.7 and a(t) ¼ 0 for other values of t. This boundary condition mimics two localized bursts of low pH, say, following two meals separated by a dimensional time interval of approximately 6 h. Results: In figure 6, we draw the evolution of the tip of level set r ¼ 0.898. For comparison, we draw this level set
J. R. Soc. Interface 11: 20140809
t
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
A three-dimensional model for the initial evolution of dental caries was considered. We took into account the diffusion of two main compounds, hydrogen ions H and calcium ions c, and the reaction at the surface of the enamel prisms. The model consists of two transport equations (2.1) and (2.2) for c and H, together with equation (2.4) for the evolution of the enamel prisms on the microscopic level.
Funding statement. This work is supported in part by an F.I.R.S.T. Initial Training Grant of the ERC, and in part by a grant from the ISF.
References 1.
2.
3.
4.
5.
6.
Fabregas R, Rubinstein J. In press. A mathematical model for the progression of dental caries. Math. Med. Biol. (doi:10.1093/imammb/ dqt008) Shellis RP, Wilson RM. 2004 Apparent solubility distributions of hydroxyapatite and enamel apatite. J. Col. Interf. 278, 325 –332. (doi:10.1016/j.jcis. 2004.06.016) Bjorndal L, Thylstrup A. 1995 A structural analysis of approximal enamel caries lesions and subjacent dentin reactions. Eur. J. Oral Sci. 103, 25 –31. (doi:10.1111/j.1600-0722.1995. tb00006.x) Bjorndal L. 2008 The caries process and its effect on the pulp: the science is changing and so is our understanding. Pediatr. Dent. 7S, S2–S5. (doi:10. 1016/j.joen.2008.02.037) Kidd EAM, Fejerskov O. 2004 What constitutes dental caries? Histopathology of carious enamel and dentin related to the action of cariogenic biofilms. J. Dent. Res. 83, C35 –C38. (doi:10.1177/154405 910408301S07) Holly FJ, Gray JA. 1968 Mechanism for incipient carious lesion growth utilizing a physical model based on diffusion concept. Arch. Oral
Biol. 13, 319 –334. (doi:10.1016/0003-9969(68) 90130-1) 7. Patel MV, Fox JL, Higuchi WI. 1987 Physical model for non-steady-state dissolution of dental enamel. J. Dent. Res. 66, 1418 –1424. (doi:10.1177/ 00220345870660090201) 8. Patel MV, Fox JL, Higuchi WI. 1987 Effect of acid type on kinetics and dissolution of dental enamel demineralization. J. Dent. Res. 66, 1425–1430. (doi:10.1177/00220345870660090301) 9. Zimmerman SO. 1966 A mathematical theory of enamel solubility and the onset of dental caries: III. Development and computer simulation of a model of caries formation. Bull. Math. Biophys. 28, 443 –464. (doi:10.1007/BF02476825) 10. Radianski RJ, Renz H. 2004 A possible interdependency between the wavy path of enamel rods, distances of Retzius lines, and mitotic activity at the cervical loop in human teeth: a hypothesis. Med. Hypotheses 62, 945–949. (doi:10.1016/j. mehy.2003.12.038) 11. Jenng YR, Lin TT, Shieh DB. 2009 Nanotribological characterization of tooth enamel rod affected by surface treatment. J. Biomech. 42, 2249– 2254. (doi:10.1016/j.jbiomech.2009.06.057)
12. Bollet-Quivogne SEC, Anderson P, Dowker SEP, Eliott JC. 2005 Microradiographical study on the influence of diffusion in the external liquid on the rate of demineralization in hydroxyapatite aggregates. Eur. J. Oral Sci. 113, 53 –59. (doi:10.1111/j.16000722.2005.00188.x) 13. Van Dijk JWE, Borggreven JMPM, Driessens FCM. 1983 Diffusion in mammalian tooth enamel in relation to the caries process. Arch. Oral Biol. 28, 391–397. (doi:10.1016/0003-9969(83)90006-7) 14. Dumont ER. 1995 Mammalian enamel prism patterns and enamel deposition rates. Scanning Microsc. 9, 429–442. 15. Hornung U. 2012 Homogenization and porous media. Berlin, Germany: Springer. 16. Bensoussan A, Lions JL, Papanicolaou GC. 2011 Asymptotic analysis for periodic structures. New York, NY: Chelsea. 17. Keller JB. 1963 Conductivity of a medium containing a dense array of perfectly conducting spheres or cylinders or nonconducting cylinders. J. Appl. Phys. 34, 991–993. (doi:10.1063/1.1729580) 18. Anderson P, Levinkind M, Elliott JC. 1998 Scanning microradiographic studies of rates of in vitro demineralization in human and bovine dental
8
J. R. Soc. Interface 11: 20140809
6. Conclusion
The equations were first solved in a rectangular geometry. Then we presented two simulations in a circular (spherical) geometry that is adequate for most initial forms of caries. We emphasized two alternative ways of quantifying caries evolution: first by visual inspection, which is modelled by following a level set of the prism density distribution, and second by measuring the total volume of enamel. The simulations performed above in circular geometry model caries progress near an occlusal surface and near an approximal surface. In both cases, the model captures the well-known triangular shapes of the global caries lesion. However, we argued that this triangular shape should be more pronounced for lesions starting at an occlusal surface. The basic model was then upgraded to incorporate the very important effect of the outer prismless layer of enamel. The extended equations model this layer as a thin domain with isotropic diffusion matrix that is resistive to dissolution relating to bulk enamel. The presence of even a very thin such layer is shown to considerably reduce the progress of caries. A very interesting feature of the model is a qualitative explanation of fig. 4 of [18]. The model presented here still contains a few undetermined parameters, such as b, R1 and g. They can be estimated by matching the model to experiments. Then the model can be used to predict outcomes of new experiments, and in general, for in silico tests of different aspects of dental caries. There are a number of ways in which the model can be further extended. One is to improve the understanding of boundary conditions. Another very interesting direction is to add into the model the effect of fluoride.
rsif.royalsocietypublishing.org
(solid line) together with the evolution of the same level set when there is no prismless outer layer (dotted line). The powerful effect of the prismless layer and its role in slowing caries evolution is clearly seen in the figure. As a further demonstration of the important role played by the prismless outer layer, we apply a boundary condition that mimics a burst of localized low pH, say following a meal, followed by a period of elevated pH that gives rise to remineralization. The geometry is rectangular as before, with the prismless enamel domain defined as x3 , j. We set equilibrium boundary condition w ¼ 0 at all boundaries except at the tooth surface x3 ¼ 0. There we set w(x1, x3 ¼ 0, t) ¼ 1024A(0.1 2 jx1j) for 0 , t , 0.2, then w(x1, x3 ¼ 0, t) ¼ 21025 for 0.2 , t , 0.4 and finally, w(x1, x3 ¼ 0, t) ¼ 0 for 0.4 , t , 0.5 to complete the cycle of approximately 8 h. Results: In figure 7, we draw the evolution of the tip of level set r ¼ 0.898 (solid line) compared with the evolution of the same level set when there is no prismless outer layer (dotted line). The effect of the prismless layer is transparent. We also note that on the scale where the graph is drawn, there seems to be a delay in the initiation of the caries evolution owing to the outer layer.
Downloaded from http://rsif.royalsocietypublishing.org/ on November 12, 2014
21. Ripa LW, Gwinnnett AJ, Buonocore MG. 1966 The “prismless” outer layer of deciduous and permanent enamel. Arch. Oral Biol. 11, 41 –48. (doi:10.1016/ 0003-9969(66)90116-6) 22. Fava M, Watanabe I, Fava de Moraes F, Da Costa L. 1997 Prismless enamel in human non-erupted decidous moalr teeth: a scanning electron microscopy study. Rev. Odontol. Univ. Sao Paulo 11, 239–243. (doi:10.1590/S0103-06631997000400003)
23. Lions PL. 1990 On the Schwartz alternating method III: a variant for non-overlapping subdomains. In Third International symposium on domain decomposition methods for partial differential equations (eds TF Chan), pp. 202 –231. Phildelphia, PA: SIAM. 24. Quarteroni A, Valli A. 1999 Domain decomposition methods for partial differential equations. Oxford, UK: Oxford Science Publications.
9
rsif.royalsocietypublishing.org
enamel. Arch. Oral Biol. 43, 649– 656. (doi:10. 1016/S0003-9969(98)00052-1) 19. Poole DFG, Shellis RP, Tyler JE. 1981 Rates of formation in vitro of dental caries-like enamel lesion in man and some non-human primates. Arch. Oral Biol. 26, 413–417. (doi:10.1016/0003-9969(81)90039-X) 20. Stumpf MPH, Porter MA. 2012 Critical truths about power laws. Science 335, 665–666. (doi:10.1126/ science.1216142)
J. R. Soc. Interface 11: 20140809
