J. MATER. CHEM., 1994, 4(5), 765-770 Statistical Analysis of Apatitic Tricalcium Phosphate Preparation Hassan Chaair, Jean-Claude Heughebaert, Monique Heughebaert and Michel Vaillant INPT, ENSCT, CNRS URA 445, Laboratoire des Materiaux, Equipe de Physico-chimie des Solides, 38 rue des 36 Ponts, 31400 Toulouse, France Apatitic tricalcium phosphate was obtained by a continuous (3-5 kg per 24 h) process using the conventional double decomposition method between an aqueous calcium nitrate solution, Ca(NO,),, and an aqueous ammonium phosphate solution, (NH,),HPO,. To check the effect of certain variables on the reaction, a fractional central composite design was set up taking six variables into account: pH, (Ca/P),,e,,,, concentration of the calcium nitrate solution ([Ca2+ I), temperature (T),duration of precipitation (R) and speed of stirring (S).The limiting factors of precipitation for apatitic tricalcium phosphate are discussed.Calcium phosphates are used as bioceramics for prosthetic application^.'-^ They are based mainly on hydroxyapatite [HAP; Cal0( PO,),(OH),] and P-tricalcium phosphate [P-TCP; Ca,(PO,),]. The difficulty with most of the conventional precipitation methods used is in obtaining well defined and reproducible solids, i.e. a solid with a given Ca/P rati~.~.~Factors governing the precipitation, pH, temperature, (Ca/P),,,,,,,, , speed of stirring, etc., are not usually precisely controlled; consequently, the solid precipitated is, in fact, a mixture of various calcium phosphates that, after heating to 900°C in air, gives a solid with approximately the desired Ca/P ratio.The purpose of this paper is to optimize the continuous synthesis of tricalcium phosphate (3-5 kg per 24 h) by looking for a possible optimum in the response surface representing the relationship between the atomic Ca/P ratio of the precipitate (ca. 1.50), and investigating the variables governing precipitation. In this work the relationship between the atomic ratio of the precipitate and six quantitative variables, i.e. pH, atomic ratio of the reagents [(Ca/P),,,,,,,], concentration of the calcium solution ([Ca" I), temperature (T), duration of precipitation (R)and speed of stirring (S),was determined by a polynomial of the second degree in a set of experiments according to a fractional central composite Experimental The precipitates were obtained by a wet process using a conventional double decomposition method between a cal- cium solution, Ca( NO,),, and an ammonium phosphate solution, (NH4)2HP04.9 A schematic diagram of the appar- atus is shown in Fig.1. A 3 dm3 reactor was maintained at constant temperature, both the reagents were introduced at the same constant flow rate and the pH was adjusted to a constant value with a pH-stat, which controlled the addition of base (or acid) as necessary. The precipitate ran out of the overflow of the reactor and was filtered, washed with de-ionized water, air-dried at 80 "C and heated to 900 "C in air. The resulting material was studied by X-ray diffraction, infrared ( IR) spectroscopy and chemical analysis.X-ray diffraction patterns were recordet at room temperature with Co-Ka, radiation (A= 1.78892 A) and a Seeman-Bohlin camera, the presence of impurities was checked by comparison with the American Society for Testing and Materials (ASTM) file." Precise Bragg angles for the samples were measured with respect to lines for NaC1, which was used as internal standard." IR spectra were recorded with a Perkin-Elmer FTIR 1600 spectrophotometer using pellets consisting of 1mg powder in 300mg KBr. The phosphorus content was t (HAP or TCP) Fig. 1 Block diagram for the continuous synthesis of calcium phosphate. R,=reactor; M =stirrer; R, = pH-stat; F = filtration system; Se =dryer; Br =grinder; and T = sieve analyzed by the colorimetric method described by Gee and Deitz.12 The relative error for phosphorus determination was ca.0.1%. Calcium was determined by atomic absorption spe~trometry,'~with a relative error of 0.3%. Statistical Analysis Experimental analysis is used in agriculture, biology and ~hemistry'~to study the empirical relationship between one or more measured responses and a number of variables. This part of the paper discusses the principles governing the construction and analysis of a fractional central composite design. Calculation: two-variable design The matrix of the complete central composite design can be written: Xl x2 -+ -+ 22 factorial design (NF)+ + -a 0 7 +a 0 I axial points 2 x 2 (N,) 0 -a 0 +a 0 0 }centre points (N,)0 0 The central composite design 22 has a total number of experiments of ~I=N,+N,+N~=~~+(~x2)+No=8+No. The theoretical model equation is therefore: YI =Po +PJl+ P2X2 +PlJ: +P22Z +P12XlX2 Note, this composition plan is not centred.It is sufficient to calculate the interactions and the quadratic variables, as with the matrix of the theoretical model: xo x1 x2 x: x; XlX2 + --+ + + + +-+ + -1+-++ + -NF+ ++ + + + + -a 0 a2 0 0 + +a 0 a2 0 0X= I+ 0-a 0 a2 0 Na+ O+a 0 a2 0 + 00 0 0 0 + 00 0 0 0 1No sum 8+No 0 0 4+2a2 4+2a2 0 The sum of X; (or of X:) equals 4+ 2a2. The coefficient parameters may be written: b =(X'X)-'X'y, and the matrix XX can be calculated: n O 0 2a2+4 2a2+4 0 2a2+4 0 0 0 0 XX= 2a2+4 0 0 0 2a2+4 4 0 SYMMETRY 2a2+4 0 4 Since X: and X2 are not centred, non-zero terms can result in X'X:2a2+4 on the first row and, for reasons of symmetry, in the first column.X: is thus not orthogonal to Xi, hence the 4 appears in the triangle at X:X; (or X;X: by symmetry). X'X is not a diagonal matrix and neither is (X'X)-'; it appears that (X'X)-' is almost diagonal since it presents numerous zero values, i.e. it is almost orthogonal. Full orthogonality is obtained leading to a diagonal (XX)-',if the variables X: and X; are replaced by centred variables. If the sum of X: is equal to 2a2+4,in n experiments, J. MATER. CHEM., 1994, VOL. 4 the new variable is named U:.The same is true for Xi: 2a2+4u;ox; -~ n Finally, a is calculated by expressing Ul orthogonally with respect to U;. The theoretical model is then written: r =Pb +Pix1+Pix2 +PilUI+Pi2G +&2XlX2 Example Let us consider a complete central composite design 22,with one experiment in the centre. In other words n=8+ No= 8 + 1=9. Therefore: ui", -2a2+4 ~ 9 u?2ux?22a2+4 -___9 hence the new variables: UfI =9Xt' -(2a2+4)=9x;' -2a2-4 The model matrix can thus be written: + --5-2a2 5-2a2 + -+ + -5-2a2 5-2a2 NF=22 -+-+ 5-2a2 5-2~~ + + + 5-2a2 5-2a2 + X= + -a 0 7a2-4 -4-2a2 0 + +a 0 7a2-4 -4-2a2 0 Na=2x2+ 0 -a -4-2a2 7a2-4 0 + 0 +a -4-2a2 7a2-4 0 + 0 0 -4-2a2 -4-2a2 0 No= 1 sum 9 0 0 0 0 0 The experiments are centred but Ut is not orthogonal to U;.By making U: orthogonal to U;: 9c (u;1)(u;2)=4(5-2a2)+4(7a2-4)(-4 -2a2) i=l +(-4 -2a2)2 = 180 -144~~-36a4 3 0 *5 -4a2-a4=0 which can be written: (a-l)(a+ 1)(a2+5)=0 where the real solutions are: a= f1, the two others being imaginary. J. MATER. CHEM., 1994, VOL. 4 Calculation with More than Two Variables On calculation, with more than two variables, orthogonality can also be written: i=l Application For k =2 (k= number of explanatory variables), and No = 1 (the same as before): (4 + 2a2)222-V2+2a2I2 4--022 + (2 x 2) + 1 =O; 9 so: (a-l)(u + l)(a + 5)2=0; real solutions: a = f1 For k = 6, and No =5 (complete central composite design): (26 + 2a2)2 (64 + 2~~)~ 26-~64-=O*a= f224+(2 x 6)+5 81 For k-p=6-1, NF=2k-p=26-1=2s=32 and N0=6 (as for this case, fractional central composite design): (26-+ 2a2)2 (32 +2~~)~ =32--O*a= +226-1 -26-' +(2 x 6) + 6 50 In the present case, a fractional central composite design is analysed in which the response (y) is the atomic Ca/P ratio of the washed solid obtained, heated to 900°C in the air [(Ca/P),,,], and the variables xj are: pH, (Ca/P),,,,,,,,, [Ca2+], T, R and S, hereafter respectively called x1, x2, x3, xq, x5, and x6.Table 1 shows the fractional central composite design presented according to the standard order; the values of the coded variable Xj are dimensionless. The values of the natural variables are summarized in Table2.The 50 experiments to be run are of orthogonal design (which means that the coefficients do not change when any model parameter changes). They are the following (Table 1). (i) The first 32 experiments belong to a 26-1 factorial fractional design; the & 1 coded values Xj were obtained by calculating: Xj =(xj-~j)/Ax The additional variable, x6, is confounded (i.e. confused) with the product XlX2X3X4Xs. (ii) The next 12 experiments were the points on the six axes, at a distance +a from the centre. (iii) In theory, we should have performed six experiments in the centre, but it is usually accepted that only three or four central experiments are sufficient to account for the coefficients of the estimated model with good accuracy (see Appendix).The distance c1 was calculated so as to obtain square values of the variables X: that are mutually orthogonal; in the present design, which consisted of six variables belonging to 26-1 and six experiments in the centre, we have a=2. The estimated model was: 6 66 6 jj =bo + C bjXj i-C 1 bjj, XjXj, + 1bjjU5 j= 1 j=1 j' = 1,j#j' j=1 Let b,X, be the general term of 9, the 28 terms (1constant + 6 variables j + 15 interactions jj' + 6 squared variables jj = 28) generally used for the construction of the model must be mutually orthogonal, and the normal equation gives the b, coefficients with the least-squares method: n b, = 7,r, where r, = c Xiuyi, i= 1c x u i= 1 Xi,and yi being the Xu and y values for the ith experiment; Y, is named contrast.Results Table 1 shows the experimental data for each atomic ratio of precipitate washed and heated at 900°C (yi). The run was performed in a randomized order (Table 1). The 28 terms were easily calculated by substituting data values in the expressions for the least-squares estimates of the coefficients (Table 3). The fitted response surface is, if expressed in real variables:-(Ca/p),,, x lo3= 1487.1+24.4Xi + *.* -7.2x6 + 2.8X1X, f ..' -3.OX5X6 + 2.9XIl-k '** + 3.7x66 (1) From this equation, it is possible to compute estimated values jji and the corresponding residuals ei = yi-9, (Table 1). An estimate of the variance of the experimental error (sf) was obtained by dividing the residual sum of squares, C,eZ (Table l), by v (number of degrees of freedom =number of experiments minus number in the model, i.e.50-28=22), see Table 4: s; =(12.44 x 10-3)/22 = 5.65 x 10-4 The statistical errors for P and Ca of 0.1 and 0.3%, respectively, agreel23l3 with a well defined product, after many measurements. But an error of s, =0.024 and relative error of 0.024/1.500 = 1.6% after only one measurement is expected (colorimetric and atomic absorption). The estimated variances of the coefficients sb2, given in Table 3 were therefore calculated by the following formulae: Ii The significance of the effects can be estimated by comparing the values of the ratio bi/si, to a critical value for the F distrib~tion,'~as indicated in Table 3.It appears that only the main effects of pH, (Ca/P),,,g,nts, [Ca"], T and S, the interactions pH-[Ca2+], T-R, (Ca/P),,,,,n,,-S, [CLi2+1-S and T-S, and the quadratic term R2 are significant. The best fitting response function is then conveniently written as follows: (cTP)whs x lo3= 1487 +24x1 + 11x2 f 7x3 + 11x4 --7x6 (f4) (f4) (+4) (f4) (f4) + 8x1 X3 +9xZx6 -SX3X6 + 7X4X6 (f4) (f4) (i4) (f4) + loX4x5 -10x55 (2)(+4) (*3) The numbers in parentheses below the coefficients represent the standard deviations (Table 3), for example: sbjj= 48.84 x = 0.0030 = 3 x lop3 Discussion Investigation of eqn. (2) showedAat if XI =0, X2=0, X3= 0, X4= 1, X, = 1 and & =0, (Ca/p),,, x 1.50 (Fig. 2). This 768 J. MATER. CHEM., 1994, VOL.4 Table 1 The fractional central composite design presented according to the standard order. (Measurements not actually carried out are in parentheses, see Appendix) order coded units of variables (Ca/p)whs exp (Ca/p)whs cal residues x3 x4 X, x6 Yi 9i lo3eiactual logical Xl x, 1 1 ------1.484 1.466 + 18 14 2 + ----+ 1.47, 1.453 + 22 37 3 -+ ---+ 1.44, 1.456 -10 21 4 + + ----1.46, 1.469 + 00 42 5 --+ --+ 1.43, 1.427 + 08 8 6 + -+ ---1.54, 1.528 + 21 36 7 -+ + ---1.46, 1.488 -20 20 8 + + + --+ 1.516 1.522 -06 29 9 ---+ -+ 1.459 1.462 -03 47 10 + --+ --1.456 1.479 -23 32 11 -+ -+ --1.45, 1.444 + 08 23 12 + + -+ -+ 1.536 1.517 + 19 13 13 --+ + --1.47, 1.457 + 18 41 14 + -+ + -+ 1.484 1.486 -02 15 15 -+ + + -+ 1.470 1.458 + 12 28 16 + + + + --1.501 1.517 -16 11 17 ----+ + 1.41, 1.403 + 16 27 18 + ---+ -1.45, 1.470 -12 35 19 -+ --+ -1.44, 1.447 + 01 18 20 + + --+ + 1.46, 1.487 -18 33 21 --+ -+ -1.440 1.459 -19 40 22 + -+ -+ + 1.44, 1.455 -09 49 23 -+ + -+ + 1.461 1.439 + 22 2 24 + + + -+ -1.550 1.548 + 02 3 25 ---+ + -1.48, 1.479 + 05 4 26 + --+ + + 1.533 1.514 + 19 16 27 -+ -+ + + 1.481 1.503 -22 48 28 + + -+ + -1.51, 1.528 -09 43 29 --+ + + + 1.42, 1.425 -01 45 30 + -+ + + -1.54, 1.539 + 09 12 31 -+ + + + -1.461 1.484 -23 25 32 + + + + + + 1.546 1.565 -19 17 33 -2 0 0 0 0 0 1.44, 1.450 -06 30 34 +2 0 0 0 0 0 1.55, 1.548 + 10 50 35 0 -2 0 0 0 0 1.43, 1.464 -34 26 36 0 +2 0 0 0 0 1.54, 1.510 + 39 38 37 0 0 -2 0 0 0 1.45, 1.463 -05 7 38 0 0 +2 0 0 0 1.500 1.490 + 10 39 39 0 0 0 -2 0 0 1.45, 1.464 -09 31 40 0 0 0 +2 0 0 1.52, 1.506 + 14 9 41 0 0 0 0 -2 0 1.41, 1.439 -24 46 42 0 0 0 0 +2 0 1.482 1.454 + 28 10 43 0 0 0 0 0 -2 1-53, 1.516 + 20 5 44 0 0 0 0 0 +2 1.47, 1.487 -15 6 45 0 0 0 0 0 0 1.49, 1.487 + 03 (44) 46 0 0 0 0 0 0 ( 1.490) 1.487 + 03 (24) 48 0 0 0 0 0 0 (1.48,) 1.487 -02 (34) 50 0 0 0 0 0 0 (1.48,) 1.487 -05 22 47 0 0 0 0 0 0 1.48, 1.487 -02 19 49 0 0 0 0 0 0 1.48, 1.487 -05 Table 2 Factorial levels used to optimize apatitic-TCP natural variables (xj) coded variables X,, X,, X,, X,, X,, Xga x1 =pH 5.25 5.75 6.25 6.75 7.25 x2 = (Ca/P)reagents 1.44 1.47 1.50 1.53 1.56 x, = [Ca2+ ]/moll-0.75 1.oo 1.25 1.50 1.75 x4 = T/"C 57.5 60.0 62.5 65.0 67.5 x5 = R/h 2.00 3.oo 4.00 5.00 6.00 x6 = S/rpm 400 500 600 700 800 OX, = (XI -6.25)/0.50; x,= (X2 -1.50)/0.03; xs = (Xj -1.25)/0.25; x4 = (X4 -62.5)/2.50 x5 = (Xg -4.00)/1.00; and x6 = (xg -600)/100.3. MATER. CHEM., 1994, VOL. 4 769 Table 3 Analysis of variable effect coefficient variance F-value (bU) (&> (b,2/&) significance 6 1.4871 ---0.0244 1.41 x 42.11 **** 0.0115 1.41 x 10-5 9.35 *** 0.0067 1.41 x 10-5 3.17 * 0.0106 1.41 x 10-5 7.95 *** 0.0036 1.41 x 10-5 0.92 NS -0.0072 1.41 x 10-5 3.66 * 0.0028 1.77 x 10-5 0.44 NS*0.0082 1.77 x 10-5 4.15 0.0037 1.77 x lo-’ 0.77 NS 0.0027 1.77 x 10-5 0.4 1 NS -O.OOO6 1.77 x 10-5 0.02 NS -0.0065 1.77 x 10-5 2.39 NS “ 3 0.0048 1.77 x 10-5 1.31 NS $ 0.0044 1.77 x 10-5 1.09 NS -0.0017 1.77 x 10-5 0.16 NS**0.0098 1.77 x 10-5 5.43 0.0022 1.77 x 10-5 0.27 NS ‘t.*0.0086 1.77 x 10-5 4.18 * 0.75 1.00 1.25 1.50 1.75-0.008 1 1.77 x 10-5 3.71 0.0074 1.77 x 10-5 3.09 * x, = [ca2+]/rnolI-’ -0.0030 1.77 x 10-5 0.51 NS Fig.2 Response function contour lines [eqn. (l)], Xl =pH =0; X2= 0.0029 8.84 x 0.95 NS (Ca/P),,,,,,, =0; X4 = T = 1; and X,= S =0. Broken line shows the o.OOO1 8.84 x 0.04 NS experimental range. -0.0025 8.84 x 0.7 1 NS -0.0004 8.84 x lop6 0.02 NS***-0.0102 8.84 x 11.77 0.0037 8.84 x lop6 1.55 NS (Fig. 2). When the calcium concentration of the C’a(N0,)2 ****: significant at a level of 0.1% (Fo.ml(1,22) = 14.38); ***: signifi-solution and the duration of the precipitation increasecant at a level of 1 % (Fo,ol(1,22) = 7.95); **: significant at a level of together, or when the duration of the precipitation increases, (1,22) =4.30); *: significant at a level of 10% (Fo.lo5% (F0.05 (1,22)= the calcium concentration of the Ca( N03)2 solution remaining 2.95); and NS: non-significant.unchanged, then the Ca/P ratio of the precipitate increases up to 1.50. Table 4 Regression variance analysis for model ratio Ca/P of the solid washed, air dried and heated to 900°C -Conclusionssource of sum of degrees of mean variation squares freedom square Fexpu Sb The precipitation of apatitic tricalcium orthophosphate, has been studied using a factorial central composite design. The regression 0.05807 27 0.002 15 3.80 -‘ response equation for the atomic ratio Ca/P of the precipitate residue 0.01244 22 0.000565 --sum 0.07051 49 ---obtained was established. From this equation, it was possible to forecast the optimal conditions to obtain 8-TCP from ‘Fexp: Snedecor factor.bS: significance test; and ‘: significant apatitic tricalcium orthophosphate, after heating to 900 “C, approximatively at a level of 0.1% (Fo,wl(27,22) = 3.85)14. and also to determine the Ca/P ratio of the precipitates. Appendix Calculations are elementary for Table 1, which show:, the six results required for the centre, this means 50 measurements. However, it is preferable to recalculate the results from the 47 measurements actually carried out. Some confounding appears only with the six square values, Xjj.So, the six normal equations have to be solved simultaneously. Because the central composite design is well balanced, the six normal diagonal coefficients of the unknown bjj are all equal to 1408/47: (-1): 1487+ ... -10(-1)’=1477 (0): g=1487 40’ (+1): Q = 1487 + ... -lo(+ 1)’ = 1477 = [( +(1)’ + ... + ( 1)2 +(4)’+ (q2]-47 1600 1408=64--=-47 47 The same holds for the 32 non-diagonal coefficients of the 7 70 unknowns bjfj,(j’# j)= -96/47: c (Xjj-XJJ)(xjrj. 47 -x,.,.)= c(XjjXjjj,)-c xjjcxjrj. j#j’ 40 x 40 1600 ~= (32 + 0) -47 = 32 --47 96-_--47 So, by addition of the six equations, we obtain a 7th equation in which all the coefficients of the left-hand member are equal to [1408 -(6 -1) x 96]/47 = 928/47.Therefore, by adding the jth normal equation to the 7th equation, with appropriate multipliers 29 and 3, the parameter bjj is isolated: bjjC(29 x 1408 + 3 x 928)]/47 =(3 + 29)xj x 3 1 Tcj, j’fj i.e. 928bjj= 32Xj + 3 1YJtj, j’#j and the corresponding contribution is: Finally, the six contributions bjj are in the j order from 1 to 6, i.e. 0.0027, -0.0001, -0.0028, -0.0006, -0.0104 and 0.0035. As in Table 3, the X,, contribution is much greater than the others. The new value of b,, is -0.0104, which is close to the -0.0102 in Table 3. Glossary b = coefficient matrix. b, = coefficient of the polynomial model. (Ca/P),,, = Ca/P ratio of washed heated solid. ei = residual of the ith experiment, ei = yi-ji.k = number of explanatory variables. p = number of complementary variables. s& =estimated variance of coefficient b,. s;? =residual variance, s;?= xie;/v. X= matrix of the model. J. MATER. CHEM., 1994, VOL. 4 xj = natural variable x for element j, and .‘cj its mean, i.e. either pH, (Ca//P)reagents, [Ca2+], temperature (T),duration of precipitation (R) or speed of stirring (S). Xj = coded variable X for element j. Xjj= centred squared variable Xj (Xjj= Xj -Xg). yi = measured response for the ith experiment. ji= calculated response for the ith experiment. Y,= contrast of a linear combination of values of y,, where the sum of the coefficient is zero. IX = distance from the centre of the design. Ax = difference between x and X.v = number of degrees of freedom = number of experiments minus number of coefficients in the model. References 1 J. G. J. Peelen, B. V. Rejda and K. de Groot, Ceramurgiu lnt., 1978, 4, 71. 2 M. Heughebaert, R. Z. Legeros, M. Gineste and G. Bonel, J. Biomed. Muter. Rex, 1988, 32(A3), 257. 3 J. C. Chae, J. P. Collier, M. B. Mayor and V. A. Surprenant, in Bioceramics: Material Characteristics versus in viro behavior, ed. P. Ducheyne and J. Lemons, Ann. New York Acad. Sci., 1988, 523, 81-90. 4 M. Jarcho, R. L. Salsbury, M. B. Thomas and R. H. Doremus, J. Muter. Sci., 1979, 14, 142. 5 J. C. Heughebaert and G. Montel, Calc. Tiss. Inr., 1982, 34, 103. 6 J. C. Heughebaert and G. Montel, Bull. Soc. C‘him. Fr., 1970, 8-9,2923. 7 G. E. P. Box and K. B. Wilson, J. R. Stat. SOC., 1951, B13, 1. 8 G. Sado and M. C. Sado, Les plans d’experiences. De l’expkrimenr- ation a l’assurance qualite, AFNOR technique, Paris, 1991. 9 J. C. Heughebaert, These de Doctorat d’Etat, INP, Toulouse, 1977. 10 ASTM, Powder diflraction jile, no. 9-1 69, fl-calcium orthophos- phate Ca, (PO4)*,ASTM, Philadelphia, PA. 11 H. E. Swanson and R. K. Fuyat, NBS, Circular 539, vol. 11, 1953. 12 A. Gee and V. R. Deitz, Anal. Chem., 1954, 25, 1320. 13 M. Pinta, Sgectrometrie d’absorption utomique, Tome 11, ed. S. A. Masson, Paris, 1972. 14 G. E. P. Box, W. G. Hunter and J. S. Hunter, Statistics for Experimenters, An introduction to Design, Data Analysis and Model Building, Wiley, New York. 1978. 15 G. E. P. Box and N. R. Draper, Empirical Model-Building and Response Surfaces, Wiley, New York, 1987. Paper 31035765; Revised 22nd June, 1993