Plaque Fluid and Diffusion: Study of the Cariogenic Challenge by Computer Modeling G.H. DIBDIN Medical Research Council Dental Group, The Dental School, Lower Maudlin St., Bristol BS1 2LY, England

Every cariogenic challenge involves a mixture of convective transport, diffusive transport, and biochemical reactions, plus physico-chemical reactions (including charge-coupling of diffusion), all of which together require numerical methods for their analysis. This presentation describes a one-dimensional finite-difference computer model of the cariogenic process, and some conclusions obtained from it. Sugar clearance from the mouth, together with site-dependent exchange between the bulk saliva and plaque surface via a salivary film, is combined with a finite-difference model of events occurring within the dental plaque. The latter includes: sugar diffusion and pH-dependent acid production; diffusion and dissociation equilibria for two acid end-products of fermentation and their anions (acetate and lactate); diffusion and dissociation equilibria of phosphate buffer; diffusion of potassium and chloride; diffusion of protons and simultaneous equilibration with fixed and mobile buffers. So that proper concentration distributions consistent with local charge neutrality can be ensured, an algorithm called Q-couple is used to impose charge-coupling between the fluxes of different ions including fixed charges. Mineral dissolution and precipitation are modeled as part of the same equilibrium calculations. The predictions of the model are compared with those of an earlier, much simpler one, in which fixed buffers were not included. It is shown that the known concentration of fixed buffer greatly extends the low pH of a Stephan curve. The isoelectric point of the plaque bacteria also appears to be of importance. The effects of various concentrations of mobile buffers, including acetate, are investigated. It is also shown that varying plaque/saliva contact over the known range derived from published studies has a profound effect on the modeled

cariogenic challenge.

My brief in this symposium is to describe computer modeling studies of those parts of the caries process which are associated with the dental plaque and plaque fluid. During a cariogenic challenge following ingestion of sugar, a complex mixture of many different processes occurs at four linked but spatially separate sites. Thus, there are simultaneous and interacting events taking place: in the bulk saliva, at the saliva/ plaque interface, in the bulk of the plaque, and in the enamel. In this presentation I am concerned mainly with events in the plaque; events in the enamel will be treated as a simple loss of mineral at the tooth/plaque interface; events in the saliva and at the saliva/plaque interface will also be considered, but only briefly. The great complexity of interactions in the plaque is such as to preclude analytical solution, but numerical methods using computer models are a useful substitute and can, I believe, give considerable insight into the processes involved. Few attempts seem to have been made in this direction. Higuchi et al. (1970) described a steady-state model which assumed constant salivary glucose and neglected intrinsic plaque buffering. More recently, a brief report by Dibdin and Reece (1984) and a full paper by Dawes and Dibdin (1986) reported time-dependent computer models, the latter including effects due to intrinsic buffering by the plaque, as well as pH-dependent acid production rates.

Materials and methods. (1) The oral cavity.-Concentration of sugar after, say, a standard 10% glucose rinse falls approximately logarithmically with time (Goulet and Brudevold, 1984; Dawes and Watanabe,

J Dent Res 69(6):1324-1331, June, 1990

(>

0

Introduction.

~C(duct) 'clearrate'

'tflmrate'

I am honored to have been asked to participate in this symposium with, among others, the two "inventors" (Edgar and Tatevossian, 1971) of the concept of plaque fluid. This concept has been, and I believe will continue to be, an extremely useful one. This notwithstanding, I want to start by emphasizing that plaque fluid can never be considered in complete isolation from the plaque as a whole. In computer modeling, one is usually concerned with breaking complex systems down into simple parts, before putting them back together again mathematically with all their interactions. Applying this procedure to plaque has been a continual reminder to me of the integral nature of this material. I hope to clarify this point of view as I proceed. Presented at a symposium on "Plaque Fluid: Biochemistry and Properties of the Plaque/Enamel Interface", held during the 18th Annual Session of the American Association for Dental Research, March 16, 1989, San Francisco, California This work was supported by the British Medical Research Council under Grant No. G8623338.

1324

C(film)

aC(f lm)

=

'fllmrate-* [ C(bulk)

-

'AC (bulk)

=

'clearrato'. [ C(duct)

- C (bulk) ]

C(fllm) ]

Fig. 1 -Model for clearance from the main salivary pool and for exchange between it and salivary films on the teeth. Concentrations in film, salivary pool, and salivary duct are assumed to be governed by the differential equations shown, where "clearrate" is a rate which defines the halving time for oral clearance, and "filmrate" is one which defines the halving time for exchange between film and salivary pool.

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

COMPUTER MODEL OF PLAQUE:CARIOGENIC CHALLENGE

Vol. 69 No. 6

1987) with a halving time of about two to three min. This is very simple to model mathematically. Concentrations of other

species will also change, with the same halving time, tending toward their duct saliva values. (2) The saliva/plaque interface.-Collins and Dawes (1987) have shown that the residual volume of fluid in the mouth would, if spread evenly over its internal surface, form a film roughly 100 Elm thick. Studies by Lecomte and Dawes (1987) and independently by Weatherell et al. (1986) both imply that this salivary film's contact with the bulk saliva is very sitedependent. This situation too can be modeled quite simply, by means of a site-dependent exchange rate between film and bulk saliva, defined with a halving time used for exchange. Fig. 1 shows the modeling scheme for both clearance and film exchange.

(3) The dental plaque.-The processes going on in the plaque itself are much more complex, and finite-difference (FD) numerical methods are appropriate. If one assumes a planar system, the imaginary plaque is divided into a series of J1 slices parallel to the tooth surface, as represented by the solid lines in Fig. 2. All chemical species and their reactions in the Jth slice are then supposed to be concentrated at its nodal plane, as represented by the broken line J, while transport is assumed to occur between neighboring nodes according to the laws of diffusion. It is possible to model a film of any chosen thickness in the scheme used, since each nodal width can be defined independently of the others. Various additional refinements are included. One of these is to define a plaque fluid volume, and make assumptions about partitioning between bacterial intracellular and extracellular (plaque fluid) compartments. In the earlier report by Dawes and Dibdin (1986), the model was simplified into (a) inward diffusion of sugar from the plaque surface; (b) microbial metabolism, converting the sugar into acetic and lactic acids (according to published apparent enzyme kinetics for plaque); and (c) diffusion and buffering of the acids thus generated. Published diffusion coefficients for acids and sugars in plaque (Dibdin, 1981, 1984) were used, which were approximately 1/3 to 1/4 of the values for aqueous solution; the buffering curves for whole plaque were those given by Stra'lfors (1950). At that time, we made no attempt to consider separate diffusion of the undissociated acid and its anion, since their diffusion coefficients are known to be very similar in aqueous solution (Albery et al., 1967); we simply solved for pH using the measured buffering of whole plaque together with the self-buffering of the computed total concentrations of acetic and lactic acid remaining at a particular point. The probable pH dependence of acid production was also included. All parameters were either taken from published figures, or else measured by us.

1325

This model predicted that simple rinsing with saturated sugar solution would give a steeply falling Stephan curve with a rather slower return to the resting pH. However, the pH minima were not, I think, as low as we had expected. Nor were the periods of pH rise as extended. Results suggested that plaque thickness would be an important variable, interacting with the acid-producing ability of the plaque. Thus, in very thick plaques, with a high metabolic activity, all the sugar appeared to be converted to acid in the outer layers, reducing the pH fall at the tooth surface. However, one of our concerns about this simplified calculation procedure was the way in which it dealt with plaque buffering. It effectively assumed all the plaque buffers to be mobile, diffusing at a rate similar to that of the undissociated acids or their anions. In fact, we know, from work on whole plaque (Shellis and Dibdin, 1988) and on single cultures (Dibdin and Shellis, 1988), that much of the plaque buffering is not easily removable. We find this to be the case even after acid washing at low pH (unpublished work). Any fixed buffering is probably due (Fig. 3) to acidic groups on the bacterial cell wall with fixed negative charges which, although they may also bind cations, are neutralized by uptake of hydrogen ions as the pH is lowered. Neutral groups which become positively charged as they take up protons are also possible. This sort of ion-exchange behavior, and indeed the presence of a multiplicity of charged species, all with the possibility of diffusing and/or interacting, has led me to develop a more detailed model for dental plaque, as described next. Such models are rather difficult to implement, however, because of strong interactions between diffusive fluxes produced by the charge-coupling. Computational methods. -The thermodynamically proper procedure for dealing with coupled ionic diffusion is to convert ionic fluxes into fluxes of neutral components (pairs of oppositely charged species). Unfortunately, this is extremely difficult to apply rigorously, and has been managed only for simple systems. A less rigorous but useful and more physically meaningful procedure is to model diffusion using the separate ion diffusion coefficients, with inter-ionic electrical forces superimposed by means of some sort of diffusion potential function. This has the additional advantage that effects of fixed charges can be modeled, an important factor when considering exchange equilibria and diffusion in the saliva/plaque/tooth system, where fixed buffers are present. The procedure described here is along these lines.

PLAQUE FIXED BUFFERING 4-

FILM

PLAQUE

- - - J = Ji

-

- --

-

,--- SALIVA POOL

$9RH

- - - J = J 1-1 - - -

R-

I

H+

R

RH

R

RH

i R

- PLAQUE - - - J = 2 ENAMEL _ - -J = 1 - - Fig. 2-One-dimensional finite-difference scheme for modeling a layer of dental plaque located on a tooth surface and covered by a salivary film. The calculation zone is divided into a series of J1 slices, with all concentrations and reactions assumed to be located at their nodal planes (dashed lines). Diffusion takes place between these nodes.

RH

H R

H~~~~~R ~ ~ H R

R

R H

Fig. 3-Schematic representation of plaque bacteria with partially protonated acidic groups R- on their outer surfaces; volume of plaque fluid bathing them greatly exaggerated.

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

J Dent Res June 1990

DIBDIN

1326

The program. -The complete program conserves both matter and charge. Flexibility is achieved by splitting up the various functions into separate routines, while efficiency of computation is enhanced by the use of a novel single-pass charge-coupling routine (Dibdin, 1989a) which I have called Q-COUPLE and have shown simulates well the effects of ionic coupling; computational efficiency is also helped by specializing as many computing functions as possible. The program is implemented on an AT-compatible microcomputer using an 80286 processor and 80287 math co-processor, in 'quickBASIC4' (Microsoft Co., Redmond, WA). Double precision was used for all floating-point computations. The program works as shown in Fig. 4. (a) Oral clearance (or sugar rinsing) plus salivalfilm exchange are separately modeled, as in Fig. 1, so that concentrations of all species in contact with the plaque are continuously updated. (b) Sugar diffusion and utilization are modeled using an iterative Crank-Nicolson central difference scheme, with timesteps of a second or less. Sugar utilization and acid production by the bacteria are computed within this iteration loop, but their dependence on pH and end-products of fermentation, being a slowly varying function, is adequately controlled by values computed at the previous time-step. Lactic and acetic acids are the only two products considered. Mobile species included (other than a sugar) are: protons, acetic acid, acetate ions, lactic acid, lactate ions, di- and mono-hydrogen phosphate, calcium, and potassium chloride. Very recently, carbonic acid and bicarbonate have been added, although the results given in this presentation do not include them. (c) Physico-chemical reactions are all assumed fast, and are computed iteratively for one node at a time, using Eqs. 10a, 10b, and 13 of the appendix. Equation 13 defines the hydrogen ion concentration in terms of dissociation constants KA, KL, KP, and total concentrations AT, LT, PT, and HT of, respecstart SUGAR RINSE

~~~has rinse

I n

Oral

finished?

o

CLEARANCE/salivary flow

EQUILIBRATE film/bulk saliva (all species)

tively, acetate, lactate, phosphate, and hydrogen in the system. The protonated fixed buffer groups are all lumped together and designated separately in Eq. 10a. This relation is approximately linear as shown, and is based on titration measurements for washed plaque samples. Equation 10b defines an iso-electric point (pH of net neutral charge) for this fixed buffer. In the concentration distributions as printed out, this is simply tabulated in terms of milli-equivalents/L of positive or negative fixed charge. The tooth surface is represented by an arbitrary (large) concentration of hydroxyapatite. Its dissolution and precipitation here and throughout the plaque are included in the equilibrium by allowing stoichiometric changes (Eqs. 23a-d) in the constituent ion totals (HT, PT, and Ca), until the solubility product is reached, with an appropriate change in mineral (HAP) concentration. These new totals, which are controlled by the solubility/dissociation equilibrium of Eq. 22 (consistent with mineral availability), are then substituted in Eq. 13. Activity coefficients have some influence on the dissociation calculations. They are computed by means of the Davies equation (Davies, 1962) using the value of ionic strength from the previous time-step. For all the ions considered, this gives activity coefficients within one to three percent of the best available values. These are incorporated into the dissociation constants and solubility products. Returning to Fig. 4, all the ions considered are then allowed to diffuse over a time-step, without any charge interactions, using the single-pass TDMA Thomas algorithm and CrankNicolson central time-difference scheme (Smith, 1985, p. 19). Acetic and lactic acids are made to appear at appropriate rates at each node where sugar has been utilized in the sugar diffusion step; the assumed ratio is 2 moles of acid/mole of glucose and 4 per mole of sucrose. For the diffusion calculations, activity corrections are neglected; when tested, they had very little effect on the results, due to the small activity coefficient gradients at the ionic strengths involved. Because the free diffusion coefficients of the ionic species differ, a small charge separation occurs over a time-step. The sonically coupled diffusion is regained, however, by following this free diffusion step with a single-pass ionic conduction step using the algorithm "Q-COUPLE" (Dibdin, 1989a) which moves them, according to their ionic mobilities and local concentrations, back to the positions they would have occupied with charge-coupling included. The whole cycle of computation is then repeated for a new time-step, after appropriate exchange of salivary film fluid with the bulk saliva, and di-

compute SUGAR DIFFUSION, with UTILIZATION of units at each node (j)

TABLE

iterate for CHEMICAL EQUILIBRIUM at each node for ions, acids etc including fixed'buffer and mineral free DIFFUSION of all ions etc. (TDMA algorithm) GENERATION of my*Qjlunits of acid at each node (j) restore charge-neutrality at each node (j) using algorithm "Q-COUPLE"

update ionic strength; activity coefficients and pH dependence of sugar utilization for each node (j)

optional PRINTOUT; STOP at t(max); MAKE TIME-STEP

Fig. 4-Flow diagram of computer modeling sugar consumed at node J re-appear as

program.

Q1 moles of

my.Q3 moles of acid

at the same

node, where my is the molar yield of acid (assumed to be 2 for glucose). Simultaneously with diffusion of sugar, the program models charge-coupled diffusion of 12 species reacting with one another, with hydroxyapatite, and with fixed charges on the bacteria.

CONCENTRATIONS AND pH ASSUMED IN BATHING SALIVA, SALIVARY FILM, AND PLAQUE AT START OF THE SIMULATION FOR FIG. 5 Dental Salivary Film Plaque 576 0 576 Sugara 7.0 7.0 7.0 pH 1 1 1, (10 or 60)b Acetatea 0 0 0 Lactatea 28.0 28.0 72.5c Potassiuma 12 12 12 Phosphatea 1.0 1.0 1.0 Calciuma 10 10 10 Chloridea 0 0 Fixed charged -34 (-68)b aConcentrations in mmol/L. The sugar concentration in the saliva persists at this level for two min before clearance commences. bNumbers in brackets show the variations for the other simulations. clnitial potassium concentration in the plaque modified where necessary to fulfill charge neutrality requirements. Bulk Saliva

dmEq/L.

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

t-. I=- AT COMPUTER MODEL OF PLAQUE:CARIOGENIC CHALLENGE

Vol. 69 No. 6

lution of this latter by mixing with duct saliva and swallowing. The printout includes a column indicating the residual charge imbalance at each node. When double-precision computation is used, this is always negligible, being less than 10-12 mol/L. The program provides screen read-out checks which are obtainable by the pressing of a pre-programed key at any time during a simulation. One of these directs the computer to a subroutine at the end of the next reaction calculation and shows all the important values for that node. Totals (HT, AT, LT, etc.) are shown, together with a comparison of assumed thermodynamic dissociation constants compared with those computed from the concentrations and activity coefficients. Tightening up the convergence criterion increases the number of iterations, and gives closer and closer match between the two. Another option allows for continuous screen read-out of the state of apatite saturation throughout the computation zone.

Results and discussion.

These are confined to studies assuming plaque thicknesses of 0.5 or 0.2 mm, and take the form of simulated Stephan curves following a 10% glucose rinse. The Table shows concentrations of the main diffusing and reacting species assumed at the start of the simulation shown in Fig. 5. The numbers in brackets relate to some of the other simulations, but it must be realized that these values also change some of the other concentrations, notably that of potassium, in order to maintain charge neutrality. Fig. 5 compares simulations using the fixed buffer model with one obtained using the old buffer model, but with identical sugar and acid diffusion coefficients, acid production rates, and buffer strengths. In both simulations, salivary film effects were neglected (i.e., perfect mixing with the bulk saliva was assumed). As can be seen, the difference is dramatic. In the mobile buffer model, the pH returns to a value of 6.8 after passing through a minimum pH of 5.5 (the lower broken line corresponds to the curve for a saturated sucrose rinse using the same model). In the fixed buffer model, on the other hand, the pH goes through a much lower minimum (4.6), and has still returned only to a value of 6.2 after two h. Nevertheless, just as in the old model, the generated acid anions (acetate and lactate) have all diffused out after a period of 25-30 min, the pH being held down by protons supplied by the fixed acidic groups. Because the free hydrogen ion activity is so low (0.01 mmol/L at pH 5), mobile buffers which ionize between the internal plaque pH and that of the external solution are needed

in order to "export" hydrogen ions and so speed up depletion of this fixed buffer. There are such mobile buffers, of course, in the form of phosphate, bicarbonate, and high-pK fatty acids such as acetate. The two models probably correspond to two extremes of what is likely to be the real situation. The advantage of the new one is that it allows for investigation of intermediate cases. The first of these mobile buffers, the H2PO4/HP04 system, is already incorporated into the model and contributes to this pH rise. Another very important mobile buffer is the bicarbonate system (Shellis and Dibdin, 1988). Although its effects were not included in the simulations described here, more recent studies (Dibdin, 1989b) suggest that its effects may be substantial, even at resting salivary levels. With its pK near that of saliva, it appears to cause a steep pH gradient near the plaque outer surface, which gradually moves back into the depths of the plaque as it depletes fixed buffers of their protons in the next layer. As pointed out by Margolis et al. (1985), a third buffer which could be important is the acetate system, which starts to operate in the lower pH range of a Stephan curve. Effect of acetate as a mobile buffer.-With its pK of 4.8, acetate could buffer substantially in the lower half of a Stephan curve. Since acetate is one of the acids considered in the model, it is simple to include various concentrations of potassium acetate in the initial plaque distribution. Fig. 6 shows the result of changing the initial acetate in the plaque from 10 to 60 mmol/L. As can be seen, the simulation does not suggest a very big difference in the depth or extent of the Stephan curve, although the simulated mineral loss is noticeably less at 60 mmol/L initial acetate. However, as modeled here, the acetate diffuses out comparatively early on in the simulation (25-30 min), so that not much is available when needed. However, we know there to be quite a lot in resting plaque, so perhaps intracellular acetate, which is also a by-product of nitrogen metabolism, stays trapped in the cells until mobilized by a pH fall. If this were so, then the anion would become available just when needed. In support of this possibility, some studies we made, for a different reason, on lactate distribution in single-strain bacterial plaques, suggest almost identical intracellular and "plaque fluid" lactate concentrations. Presumably it is only the pH which differs in the two compartments. This is another situation where it may be essential to consider plaque fluid together with the rest of the plaque. Effect of plaque thickness.-Fig. 7 shows more Stephan 7

7

J

1327

pH

I=

6

at

pH at enamel surface

enamel surface

.4

5

.4

-

1.'.-'

0

f.-.410

==I_

4

§T

*

4

5

4 20

30

.g-.

40

50

-

-m-..

60

70

80

90

0

.

time time 100 110 120 in minutes

eqf

-.

Fig. 5-The solid line shows a Stephan curve computed using the new fixed buffer model (initial values as given in the Table) after a two-minute 10% glucose rinse. This compares with the upper dashed line which was obtained using the model of Dawes and Dibdin (1986) for, where applicable, identical conditions. The lower dashed line reproduces that obtained in Fig. 1 from the same paper, following a saturated (3.3 mol/L) sucrose rinse. Plaque thickness = 0.5 mm throughout.

10

20

30

40

50

60

70

80

90

| 100 110 120

0

mineral losns

5

time in -minutes i l- Rae

EE

at

1 surfaceUEh

enamel

(arb. units)

-

15

Fig. 6-Effect of initial acetate concentration on Stephan curve (upper graph) and simulated mineral loss (lower graph). Plaque thickness = 0.5 mm; "filmtime" = one min.

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

, ; -4

1328

DIBDIN

J Dent Res June 1990

curves, this time for a 0.2-mm-thick plaque. As expected, the pH changes are much more rapid; such rates depend inversely on the square of the thickness, and so should be about six times faster than for a 0.5-mm-thick plaque. Effect offixed buffer strength (BF). -This is also dealt with in Fig. 7. Clearly, the cariostatic role of fixed buffers is not so direct as with mobile buffers such as acetate, phosphate, or carbonate. Although these buffers will initially slow down a pH fall by taking up protons, once protonated they will act as a fixed reservoir which retards a subsequent pH rise. Fig. 7 shows this for the case of two values of the fixed buffer slope, namely, 34 and 51 Wxmol/g/pH unit. As can be seen, the pH minimum is lower and is reached more rapidly at the lower value. On the other hand, the rate of return toward neutral pH is also faster. Under these conditions, the net result (Fig. 7) seems to be that there is less mineral loss with the higher buffering capacity, although I don't think that this would necessarily always be the case. As in Figs. 5 and 6, the iso-electric point assumed here was 6. Effect of iso-electric point (iso-pH). -This is one variable for which I have, as yet, no information whatever. I have therefore chosen values of either 6 or 5 quite arbitrarily, for the purposes of modeling, and looked at their effect on (a) the Stephan curve and (b) the distribution of potassium. Fig. 8 shows the modeled Stephan curve. It will be noticed that the

7

pH at enamel surface

6

TV +

4

time

0

10

30 -40

20

50

60

70

80

90

100

110 120

--==-==--

in inute s

m

mineral loss

at enamel surface (arb. units)

10

--;;Jim 1 & i;

1

=

=

m

=~~~-.-

L

=

pH profile goes to a slightly lower minimum with the lower iso-pH, and the subsequent pH rise is slower and approaches a lower asymptotic value. Another consequence of charge-coupling between diffusive fluxes of ions is that certain of them may be forced to move against their own concentration gradients. The lower curve was produced assuming a salivary potassium content of 28 mmol/ L, as denoted by the horizontal line. The two curves show how the potassium content might change in the inner part of a plaque during the course of the Stephan curve shown above. An isopH of 6 leads to an equilibrium potassium content in the inner plaque of 40 mmol/L at the end of a Stephan curve. Changing the iso-pH to 5 results in an increase in this equilibrium potassium content to about 58 mmol/L. It is tempting to relate this to measured potassium levels in plaque fluid. However, the plaque fluid removed by centrifugation of whole plaque would not be expected to carry with it any cations needed to neutralize fixed charges remaining on the bacterial residue. In the program as presently constituted, no provision is made for exchange of calcium with protons at the fixed charge sites. However, my colleague Keith Rose is studying calcium affinities of bacteria and plaque. It will be interesting to include such effects in the future. At the moment the simulation assumes that there is insufficient calcium to compete for sites to any great extent with the much higher concentrations of monovalent cations such as potassium. Last, I want to go back to where I started, and mention briefly the effect of the salivary film. In Fig. 9 the top graph shows a typical family of Stephan curves modeled for a 0.2mm-thick plaque after a two-minute 10% glucose rinse. The three curves correspond to three situations. The top curve corresponds to perfect contact between the salivary film and bulk saliva. The value of "filmtime" for the lowest curve was obtained by analysis of the data of Lecomte and Dawes (1987), and corresponds to the situation at a buccal site where contact with the bulk saliva is poorest. The middle curve corresponds to an intermediate situation. As can be seen, the degree of contact appears from the simulation to have a profound effect on the shape and extent of a Stephan curve. The bottom Fig. models predicted dissolution at the tooth surface for the same three situations. A point I want to stress is that most of the important variables which I have used in this modeling program have actually

if Us

Fig. 7-Effect of a change in buffer slope (fixed buffer model; BF in mEq/g/pH) on Stephan curve (upper graph) and mineral loss (lower graph). Thinner plaque results in a more rapid change.

:-A.

_t__,

4

___

pH at enamel

surface

pH

5~

at enamel surface

0

|.

i.

10

0

20

'

iI.

30

__I

40

50 60

._-_--

0

time

90

--

100 110 120 I

i-

in

minutes

potassium

conc.

in plaque

(mmol/L)

+minutes

, ti mineral 64

I-r-

70 80

4 ~~~~~~~~~~~~Itime 10 20 30 40 50 60 70 80 90 100 110 120 in

40

-=

20-

0 Fig. 8-Effect of iso-electric point of fixed buffers on Stephan curve (upper graph) and potassium balance (lower graph) between plaque fluid and saliva.

oss

at enamel surface

.m

----i =

20

,

I

10 4-

e

--

----

--- -

-

+-===:=T-

I

^

_

.-.

_-

a__=

= =

_

(arb. units) 40+

-

-

Fig. 9-Effect of halving time for exchange ("filmtime") between salivary film and bulk saliva on the shape of a Stephan curve (upper graph), and simulated mineral loss (lower graph). A "filmtime" of 4 corresponds to conditions at a fairly stagnant buccal site, while one of 0 corresponds to a site in direct contact with a well-stirred pool of saliva.

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

COMPUTER MODEL OF PLAQUE:CARIOGENIC CHALLENGE

Vol. 69 No. 6

been measured in one way or another on dental plaque, either by me or by others. An important exception to this at the present time is the iso-electric point we have assumed for dental plaque. We intend to look into this in the near future using concentration cell and other methods. There are inevitably many approximations and omissions in the model which I have described, and I've mentioned some of them. However, I would like to mention one or two others. One concerns the initial conditions assumed in these simulations (see Table). Since producing the simulations reported here, I have investigated the effect of allowing a period of equilibration between plaque and saliva before starting a simulated glucose rinse. As a result of Donnan membrane effects (due to the assumed plaque fixed charge), this results in a reduced initial potassium content in the plaque, compensated for by a corresponding reduction in mobile anions such as phosphate and chloride. However, the effect on the modeled cariogenic challenge is quite small, since most of the redistribution anyway occurs in the first few minutes of the cariogenic challenge. Nevertheless, this is an area which needs further investigation. Another important approximation concerns the assumption of instantaneous equilibrium with a single (hydroxyapatite) mineral phase, since it is known that direct hydroxyapatite deposition occurs only with difficulty, if at all. Because of this, I plan to include modeling of a precursormineral phase, probably dicalcium phosphate dehydrate. However, it is unlikely that general conclusions regarding the cariogenic challenge will be greatly affected, since buffering due to the mineral itself is usually small in comparison with that due to the other buffers. The bicarbonate buffer system, which was omitted from these simulations, has now been included in the model, and a report of its effect has recently been given (Dibdin, 1989b, 1990). Finally, studies which we have in progress on competition between calcium and protons for the fixed buffer groups should allow for more realistic modeling of the effect of calcium binding. A factor which has not been addressed at all as yet is the effect of possible intrinsic pHrise factors on the progress of a cariogenic challenge. Despite these caveats (and in a sense because of them), I hope I have convinced you that even imperfect models of the sort I have described can help toward a greater insight into the factors involved in a process such as dental caries.

Appendix: Fast reaction processes with and without mineral dissolution. The complete program works by modeling diffusion steps including acid production alternately with chemical reaction, changes due to one being held at zero while those due to the other are computed. During a reaction step, each node is therefore treated as an isolated system, the only change in total concentration of a solute being due to dissolution or precipitation of mineral. It is assumed that these reactions are fast in relation to diffusion processes, and are therefore in equilibrium. Formation of complex ions is neglected. Except to define the ion product of hydroxyapatite, phosphate is assumed to exist entirely as H2PO4- (written "H2P") and HP042- ("HP2-"). The following were the thermodynamic dissociation constants assumed for, respectively, acetate, lactate, and phosphate (2nd dissociation constant). They and the other constants quoted (Eqs. 19a-c) are for 370C and were collected by Shellis (1988) from various sources. TKA = 1.73-10-5 (la) TKL = 1.32-10-4 (lb) TKP 6.59 10-8 (lc) =

-

1329

Considering acetic acid as an example, its thermodynamic dissociation constant is: TKA

(A-)-(H+)

=

(2)

HA Round brackets here refer to activities and square brackets to concentrations. In terms of concentrations, and making the usual assumption that the activity coefficient of all uncharged species is 1, we have:

f12[A-]-[H+] [HA]

TKA

(3)

Here f1 is the monovalent activity coefficient, which we calculate from the Davies extension of the Debye-Hickel equation

(Davies, 1962): log1o(f

lO5)

[(1

-

0.31

(4)

where Z is the valency and I the ionic strength. Defining an operational dissociation constant (KA) for acetic acid in terms of concentrations, we have: KA

=

[A-J-[H+] [HA]

(5)

TKA/fl2

(6)

giving: KA

=

Re-arranging 5 and adding 1 to both sides gives: KA + [H+] [A-] + [HA]

[H+]

[HA]

Since, neglecting ion-pair formation, the numerator on the righthand side can be equated to AT, the total concentration of acetate, invariant during the reaction step, we have: KA [HH] +

[HA] = AT

(7a)

Similarly, re-arranging appropriate dissociation equations for the other acids in terms of the total amounts of lactate (LT) and phosphate (PT) present, we get:

[H+] LT [HL] [HL] = LTKL + [H+]

(7b)

KP + [H+]

(7c)

=

[H2P-] = PT [HP2-] = PT.

where it can

K + [H+]

similarly be shown that: KL TKL/f12 =

KP In

=

TKP2/f14

(7d) (8a) (8b)

addition, fixed buffering in the plaque:

H+ + R- =HR (9) associated with the bacterial surfaces, found to be fairly linear with pH, can be approximated as: HR = BF-(A - pH) (1Oa)

QR = BF-(pHi,0 - pH)

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

(lOb)

1330

DIBDIN

J Dent Res June 1990

where A is a constant, chosen arbitrarily to be greater than the highest pH envisaged; BF is the buffer slope, obtained from practical measurements of washed plaque or bacterial suspensions (pLmol/g/pH unit); pHi10 is the mean iso-electric point of the plaque; and QR is the pH-dependent fixed charge

(mEq/L). The total amount of

hydrogen (HT) in this system, including hydrogen bound to these fixed buffers (R-) to give HR, is here assumed constant. We then have: [HA] + [HL] + 2- [H2P-] + [H+] + HR = HT - [HP2-] (11)

Substituting from (Eqs. 7a-d) into (11) gives: [H ]*

2-PT

+KP + [H+] HT

=

-

HR

[H+]

[HP2-] PT-KP KP+ [H+]

[H+] =

AT

KA+ [H+]

+

-

LT

KL+ [H+]

(12) 1b (12b)

(13) 2-PT HR +1 + KP + [H+] [H+]

This value of [H+] is iterated to convergence, the values on the righthand side of Eq. 13 being replaced with the new estimate obtained from the lefthand side after the previous iteration. Once Eq. 13 is converged, the appropriate equation (Eqs. 7a-e) is used to calculate the concentration of each undissociated acid, from which each anion concentration is obtained by subtraction from the total, in order to conserve matter. Inclusion of apatite dissolution and precipitation. -When this is to be included, it is necessary to consider the ion product at each node, and to allow for corresponding stoichiometric increases or decreases of HT, PT, and [Ca2+] in solution as the mineral enters or leaves solution. These changes are incorporated within the iteration loop in Eq. 13 as follows. The hydroxyapatite dissolution equation is:

(14)

(15)

but

TKP3-(HPO2-)

4(Hf) (HPO-) =TKP2- (H2PO4-)

(OH-)-(H+)

(16) (6 (17)

~~(Hf)

4

=

Re-writing each activity in terms of concentration and activity coefficient: (20a) (Ca2+) = [Ca2+] . f14 (20b) (H2P04-) = [H2PO4-] * f (20c) (H+) = [H+] . f then:

TKW

[H2P04. TKP33 TKP23 * TKW * f,'6 DS _ [

+[H~]

Ca5 (OH)(PO4)3 + H+ = 5Ca++ + 3PO3- + H20 The thermodynamic ion activity product (IAP) is: IAP = (Ca2+)5 . (PO3-)3 - (OH-)

(19)

(21)

+ 1]

PT-KP KP +

5. (H2PO43 .TKP33 TKP23 TKW

We define the degree of saturation as:

giving: HT

P (Ca2=

1AP=

+[H+] + LT KA +H KL+ [H+]

= HT -

Substituting from Eqs. 16-18 in Eq. 15 gives:

(18)

where TKP2 (already defined) and TKP3 (6.92- 10-13) are, respectively, the second and third thermodynamic dissociation constants of phosphoric acid, and TKW (2.49-10-14) is the thermodynamic dissociation constant of water at 37°C.

AP

(22)

where TK(HAP) is the solubility product of hydroxyapatite, assumed to be 7.36- 10-60, and n is the number of constituent ions which go to make up the stoichiometric formula of the solid (9 in the case of hydroxyapatite). If Q moles/L of HAP dissolve during a time-step, then stoichiometry (Eq. 14) dictates that: PT = PTO + 3Q (23a)

(23b) (23c) (23d) [HAP] = [HAP]o- Q where the suffix o (old) indicates concentrations before any dissolution or precipitation has been allowed. We now allow Q to vary during each iteration until, as well as Eq. 13 being satisfied, either DS (Eq. 22) is unity or no solid mineral is left. Noticing (Eqs. 21-22) that DS varies with [Ca2+]519, we allow Q to change by an amount: AQ = 0.2-[Ca2+] (DS- 1-5 - 1) (24) where the values on the righthand side of the equality are from the previous iteration. This gives good convergence, and assuming that it is reached after n iterations, we get: (25) Q = AQ1 + AQ2 +.AQ. This Q is now used to calculate the new values in Eq. 23 (a-d). Last, after leaving the iteration loop, but before moving on to compute reaction at the next node, one more calculation is necessary. Although all other totals such as AT, LT, and PT are conserved to the full computational accuracy of the computer, the total value of HT at this point is not. This is because Eq. 13 is satisfied only to within the convergence criterion. For this to be rectified, the dissociation balance of one acid/ base pair (acetic acid/acetate, say) is adjusted (see, e.g., Eq. 6), until HT has a value which conserves total hydrogen within the complete solid/solution system. The consequence of this is a slight shift in that particular dissociation balance which shows up as a small error in the dissociation constant obtained by back-calculation from the final concentrations. Since this shift does not lead to cumulative errors, it is far less important than cumulative loss of material. Furthermore, its magnitude can be reduced to any required level simply by tightening up the convergence criterion. [Ca2+] = [Ca2+]o + 5Q HT HTO - Q

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.

Vol. 69 No. 6

COMPUTER MODEL OF PLAQUE:CARIOGENIC CHALLENGE

Acknowledgments. I thank Drs. Colin Dawes and Peter Shellis for many useful

discussions.

REFERENCES ALBERY, W.J.; GREENWOOD, A.R.; and KIBBLE, R.F. (1967): Diffusion Coefficients of Carboxylic Acids, Trans Faraday Soc 63:360-368. COLLINS, L.M.C. and DAWES, C. (1987): The Surface Area of the Adult Human Mouth and Thickness of the Salivary Film Covering the Teeth and Oral Mucosa, J Dent Res 66:1300-1302. DAVIES, C.W. (1962): Ion Association, 2nd ed., London: Butterworths. DAWES, C. and DIBDIN, G.H. (1986): A Theoretical Analysis of the Effect of Plaque Thickness and Initial Sucrose Concentration on Diffusion of Sucrose into Dental Plaque and its Conversion to Acid during Salivary Clearance, J Dent Res 65:89-94. DAWES, C. and WATANABE, S. (1987): The Effect of Taste Adaptation on Salivary Flow Rate and Salivary Sugar Clearance, J Dent Res 66:740-744. DIBDIN, G.H. (1981): Diffusion of Sugars and Carboxylic Acids Through Human Dental Plaque in vitro, Arch Oral Biol 26:515523. DIBDIN, G.H. (1984): A Brief Survey of Recent in vitro Work on Diffusion of Small Ions and Molecules in Dental Plaque. In: Cariology Today, Int. Congr., Zurich (1983), B. Guggenheim, Ed., Basel: Karger, pp. 191-198. DIBDIN, G.H. (1989a): A Simple Add-on Algorithm to Extend Onedimensional Finite Difference Diffusion Calculations to Include Charge-coupling, Comput Appl Biosci 5:19-26. DIBDIN, G.H. (1989b): Computer Modelling of Stephan Curves with Saliva/Plaque Exchange via a Thin Salivary Film, Caries Res 23:441, Abst. No. 55. DIBDIN, G.H. (1990): Effect on a Cariogenic Challenge of Saliva/

1331

Plaque Exchange via a Thin Salivary Film Studied by Mathematical Modelling, Caries Res (in press). DIBDIN, G.H. and REECE, G.L. (1984): Computer Simulation of Diffusion with Reaction in Dental Plaque, Caries Res 18:191-192, Abst. No. 110. DIBDIN, G.H. and SHELLIS, R.P. (1988): Physical and Biochemical Studies of Streptococcus mutans Sediments Suggest New Factors Linking the Cariogenicity of Plaque with its Extracellular Polysaccharide Content, J Dent Res 67:890-895. EDGAR, W.M. and TATEVOSSIAN, A. (1971): The Aqueous Phase of Plaque. In: Tooth Enamel U, R.W. Fearnhead and M.V. Stack, Eds., Bristol, UK: John Wright, pp. 229-232. GOULET, D. and BRUDEVOLD, F. (1984): Salivary Glucose Clearance after Rinsing with Solutions of Different Concentrations of Glucose, Caries Res 18:481-487. HIGUCHI, W.I.; YOUNG, F.; LASTRA, J.L.; and KOULOURIDES, T. (1970): Physical Model for Plaque Action in the ToothPlaque-Saliva System, J Dent Res 49:47-60. LECOMTE, P. and DAWES, C. (1987): The Influence of Salivary Flow Rate on Diffusion of Potassium Chloride from Artificial Plaque at Different Sites in the Mouth, J Dent Res 66:1614-1618. MARGOLIS, H.C.; MORENO, E.C.; and MURPHY, B.J. (1985): Importance of High pKa Acids in Cariogenic Potential of Plaque, J Dent Res 64:786-792. SHELLIS, R.P. (1988): A Microcomputer Program to Evaluate the Saturation of Complex Solutions With Respect to Biominerals, Compute Appl Biosci 4:373-379. SHELLIS, R.P. and DIBDIN, G.H. (1988): Analysis of the Buffering Systems in Human Dental Plaque, J Dent Res 67:436-446. SMITH, G.D. (1985): Numerical Solutions of Partial Differential Equations. Finite Difference Methods, 3rd ed., Oxford: Clarendon Press. STRALFORS, A. (1950): The Bacterial Chemistry of Dental Plaques, Odont Tidskr 58:155-341. WEATHERELL, J.A.; STRONG, M.; ROBINSON, C.; and RALPH, J.P. (1986): Fluoride Distribution in the Mouth after Fluoride Rinsing, Caries Res 20:111-119.

Downloaded from jdr.sagepub.com at WAYNE STATE UNIVERSITY on March 27, 2015 For personal use only. No other uses without permission.