## Abstract

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 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 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 1*a*. The outer layer of the tooth consists of enamel in the form of a packed array of rods (prisms) with radius of approximately 2−3 µm. The prisms vary considerably in length (20−1000 µm, [10]), and the inter-prism distance is approximately 1−3 µm [10,11]. The prism long axis is in general orthogonal to the tooth surface. A sketch of the prisms is depicted in figure 1*b*.

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 µm, 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 and , respectively, where and denote length and time, respectively, the term has the dimension of length. To understand the physical meaning of *λ*, consider a simple one-dimensional model for diffusion and reaction in equilibrium, given by the equationwhere *u* is the concentration of some substance. The solution is *u*_{0}e^{−x/λ}. Thus, *λ* represents the size of the interval that is most affected by the variation of the boundary concentration.

Using *D* = 5 × 10^{2} μm^{2} s^{−1} and *R* = 5 × 10^{−2} l s^{−1} [12,13], we obtain *λ* ∼ 100 µm. The small ratio l/*λ* implies a separation of scales. In what follows we shall use a non-dimensional formulation, where length is scaled by *λ* 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*(*x*_{1}, *x*_{2}, *x*_{3}, *t*) and the hydrogen ion concentration *H*(*x*_{1}, *x*_{2}, *x*_{3}, *t*). In addition, we assume that the prisms essentially retain their cylindrical shape [14]; thus, the model includes an equation for the function *ρ*(*x*_{1}, *x*_{2}, *x*_{3}, *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 inter-prism region. We next used the separation of scales, i.e. 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*/*λ*. 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 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 *x*_{3} 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 2.1 2.2

Here, and *∂ _{t}* denote partial differentiation, and

*d*is the ratio of the diffusion coefficients of the hydrogen and calcium ions. The reaction term is proportional to

_{H}*c*−

_{s}*c*. The equilibrium concentration

*c*is modelled by 2.3

_{s}The first term on the right-hand side in equations (2.1) and (2.2) models the diffusion along the ‘easy’ direction *x*_{3}. The second term models the diffusion in the direction orthogonal to the prism long axis. The effective diffusion tensor can be 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 |*Ω*|, which is the local volume per unit length of the inter-prism domain, and |∂*E _{p}*|, 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

*ρ*(

*x*

_{1},

*x*

_{2},

*x*

_{3},

*t*), the local prism radius, which is also scaled by

*l*, such that the maximal radius is 1. Note that when

*ρ*→ 1, the lateral diffusion approaches zero. The geometric parameters in this setting are

*|*

*Ω*

*|*= 4 −

*π*

*ρ*

^{2}and

*|∂E*= 2

_{p}|*π*

*ρ*. The function

*ρ*itself evolves according to the dissolution process 2.4

The model still contains three parameters: *R*_{1} is determined by the local reaction rate, *γ* is determined by the precise model for the equilibrium concentration *c _{s}* and

*b*is related to the molar concentration of the hydroxyapatite crystals. However, we prefer to leave

*R*

_{1},

*γ*and

*b*as free parameters that can later on be adjusted either from theoretical considerations or by fitting the model to observations.

As is readily seen from equations (2.1), the model is anisotropic, as the diffusion coefficient in the *x*_{3} direction (scaled to be 1 in our non-dimensional set-up) is different from the diffusion coefficients in the (*x*_{1}, *x*_{2}) directions, given by the effective diffusion matrix . This effective tensor can be computed separately from the transport equation on the microscopic scale [1]. As the microscopic geometry is fully encoded in the function *ρ*, this can be done in principle. However, practically, the problem of evaluating at any given point (*x*_{1}, *x*_{2}, *x*_{3}, *t*) is a formidable task. Fortunately, we found a way to greatly simplify this task; we computed for many representative values of *ρ*, and found out that, except for the limit of highly packed prisms, can be approximated by a linear function of , where *I _{d}* is the 2 × 2 identity matrix. We note that this simple formula does not hold when

*ρ*is near 1, i.e. at highly packed prism configuration where the prisms nearly touch each other. At extreme values of

*ρ*, a separate calculation should be done, solving in full the local equations that govern [1]. Alternatively, asymptotic formula as those derived by Keller [17] can be used.

We point out that modelling the progress of caries by following *c*, *H* and *ρ* 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 *ρ*.

## 3. Simulations in a three-dimensional rectangular geometry

In this section, we consider a rectangular tooth sample. The domain is

Before presenting specific simulations, we point out that the model above can be simplified if we approximate the factor *d _{H}* by 1. Under this approximation, we can define the function . The model then becomes
3.1together with the equation for the evolution of

*ρ*3.2We 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

*d*∼ 1 is valuable.

_{H}We assume the boundary condition *w* = 0 at all the boundaries, except for the square which is the tooth outer boundary, where we apply the boundary condition *w* = 10^{−4}. 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: *R*_{1} = 1, *γ* = 1 and *b* = 3. The initial conditions are *w* = 0 and *ρ* = 0.9.

*Results*: We present the progression of caries in two ways. First, we follow the tip of a certain level set of *ρ*. 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 *x*_{3}(*t*) of the level set *ρ* = 0.86 is depicted in figure 2*a*. Recall that the basic time unit is 20 s, so the curve follows a period of about 10 h. The total amount of enamel, determined by the total volume is depicted in figure 2*c*. 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 proportional to , where *d* is a diffusion coefficient, then caries should propagate at rate proportional to . Therefore, we depict in figure 2*b* the motion of the tip *x*_{3} as a function of . The curve seems not far from linear, but as we argue below this shape is misleading.

## 4. Simulations in a circular geometry

To simulate the model in a more realistic geometry, as sketched in figure 1*a*, 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 1*a*). In both cases, we can work locally in a spherical geometry (*r*, *θ*, *ϕ*). 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
4.1where again we use the simplified model with *w* = *H* − *c*/*γ*.

As a further simplification, and also as another example of the model in a curvilinear geometry, we interpret the drawing in figure 1*a* as a two-dimensional sketch and solve the equation in a polar coordinate framework (*r*, *θ*). The equation for *w* takes the form
4.2We solve equation (4.2), together with the equation for *ρ* and using the formula for |*Ω*| and |∂*E _{p}*| in several set-ups. In all of them, the geometry is an angular sector
4.3We simulate caries progress when it starts at an

*occlusal surface*. The boundary conditions are

*w*= 0 at all the boundaries, except at the sector

*Γ*= {(

_{s}*r*,

*θ*)|

*r*=

*r*

_{1}, |

*θ|*<

*θ*

_{0}}, where

*w*(

*r*

_{1},

*θ*,

*t*) = 10

^{−4}. We use the parameters

*r*

_{1}= 1,

*r*

_{2}= 3,

*β*=

*π*/4,

*θ*

_{0}=

*π*/30,

*γ*= 1,

*b*= 3 and

*R*

_{1}= 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 *x*_{3}(*t*) of the tip of a level set of *ρ*, as an indication of a visible lesion (figure 3*a*), as well as the total amount of enamel (figure 3*d*). The tip's evolution is also drawn as a function of to examine the importance of the diffusion process in the caries progress. We also draw in figure 3*c* the boundary of a given level set *θ*(*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 3*c* 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 of . 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 favour of the 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 2*b* and 3*b* are not convincing that the lesion propagates at a rate proportional to . 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 2*c* and 3*d*).

To demonstrate the danger in such power laws, we simulated a case where the evolution of the tip of a level set of *ρ* can be solved analytically. We use the same model and geometry as in the last example, except that the boundary condition on the curve *Γ*_{1} (defined by *r* = *r*_{1}) is now *w*(*r*_{1}, *θ*, *t*) = 10^{−4} cos(2*θ*). In light of the relatively slow enamel dissolution process, we solve the equation for *w* at a fixed geometry where *ρ* = 0.98. At this extreme value of *ρ*, we cannot use anymore the linear formula for . Instead, we evaluated the effective diffusion coefficient by solving the microscopic cell problem introduced in [1] and found . In addition, we used the parameters *R*_{1} = 1/3, *γ* = 3 and *b* = 1/3. For this choice of parameters, the equations for *w* and *ρ* are approximately
4.4The function *w* reaches a steady state that can be found explicitly
4.5where *K*_{1} is the modified Bessel functions of the second kind. It is convenient to define *ρ* = *ρ*_{0} + 10^{−4}*ρ*_{1} where for the parameters we used *∂*_{t}*ρ*_{1} = −10^{4}*w*. Substituting the expression given by *w* in equation (4.5) into the equation *∂*_{t}*ρ*_{1} = −10^{4}*w* and solving for *ρ*_{1} at the tip, that is for *θ* = 0, we obtain
4.6It follows that the tip of the level set *ρ*_{1} = 10 is given by the relation *t*(*r*) = (*K*_{1}(2)*/K*_{1}(2*r*))10.

*Results:* We use the formula for *t*(*r*) to draw the theoretical inverse curve *r*(*t*) in figure 4*c*. This curve can now be compared to the numerical solution *r*(*t*) given in figure 4*a* and the curve depicted in figure 4*b*. It is clear from these figures that the is not a correct power law.

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 *Γ*_{2} defined by *Γ _{s}* = {

*r*,

*θ*)|

*r*=

*r*

_{2}, |

*θ*| <

*θ*

_{0}}, where

*w*(

*r*

_{2},

*θ*,

*t*) = 10

^{−4}. We use the parameters

*r*

_{1}= 1,

*r*

_{2}= 3,

*β*=

*π*/4,

*θ*

_{0}=

*π*/30,

*bγ*= 3 and

*R*

_{1}= 1.

*Results:* As in the previous examples, we follow the caries progress by drawing (figure 5*a*) the function *x*_{3}(*t*), the tip of the level set *ρ* = 0.86 and (figure 5*c*) the function which is the total volume of enamel. In figure 5*b*, 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 3*c*).

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 . We follow the same level set and depict its width in figure 5*d*. Comparing figure 5*b* with 5*d* 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
5.1

The thickness of the outer layer, that we denote here *ξ*, is reported to be in the range of 15–200 µm [6,21,22]. For the purpose of our simulation below, we take the layer to be of size 70 µm, which is approximately 0.7 in the non-dimensional length unit scaled by *λ*.

We could not find in the literature direct estimates for the diffusion coefficient *D _{r}* in the prismless layer. However, the experiments of Anderson

*et al*. [18] give some indication for the magnitude of

*D*. 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

_{r}*D*∼

_{r}*ξ*

^{2}/

*T*, and taking

*T*= 10 h and

*ξ*= 70 µm, we obtain

*D*∼ 1.5 × 10

_{r}^{−1}µm

^{2}s

^{−1}. 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 non-dimensional diffusion coefficient by

*d*, which is the ratio between

_{r}*D*and the main diffusion coefficient in bulk enamel that we took earlier to be approximately 5 × 10

_{r}^{2}µm

^{2}s

^{−1}; we thus obtain

*d*∼ 1/3000.

_{r}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 in all the equations through *t* = *t*′/*d _{r}*. 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:
5.2and
5.3

Note that in the prismless layer, the geometry, i.e. *Ω*, 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 *x*_{3} = *ξ*. In addition to equations (5.2) and (5.3), we need to follow the evolution of *ρ*, although only in the bulk enamel *x*_{3} > *ξ*. In the new time scale, this equation reads
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 *x*_{3} = *ξ*. 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 *x*_{3} < *ξ*. We set as before equilibrium boundary condition *w* = 0 at all boundaries except at the tooth surface *x*_{3} = 0. There we set
5.5where *A* is the Heaviside step function and *α*(*t*) = 1 for 0 < *t* < 0.1, 0.5 < *t* < 0.7 and *α*(*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 *ρ* = 0.898. For comparison, we draw this level set (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 *x*_{3} < *ξ*. We set equilibrium boundary condition *w* = 0 at all boundaries except at the tooth surface *x*_{3} = 0. There we set *w*(*x*_{1}, *x*_{3} = 0, *t*) = 10^{−}^{4}*A*(0.1 − *|x*_{1}*|*) for 0 < *t* < 0.2, then *w*(*x*_{1}, *x*_{3} = 0, *t*) = −10^{−}^{5} for 0.2 < *t* < 0.4 and finally, *w*(*x*_{1}, *x*_{3} = 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 *ρ* = 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.

## 6. Conclusion

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.

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*, *R*_{1} and *γ*. 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.

## 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.

- Received July 22, 2014.
- Accepted August 27, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.