Effect of Process Variables on the Drucker–Prager Cap Model and Residual Stress Distribution of Tablets Estimated by the Finite Element Method

Abstract

A multivariate statistical technique was applied to clarify the causal correlation between variables in the manufacturing process and the residual stress distribution of tablets. Theophylline tablets were prepared according to a Box–Behnken design using the wet granulation method. Water amounts (X₁), kneading time (X₂), lubricant-mixing time (X₃), and compression force (X₄) were selected as design variables. The Drucker–Prager cap (DPC) model was selected as the method for modeling the mechanical behavior of pharmaceutical powders. Simulation parameters, such as Young’s modulus, Poisson rate, internal friction angle, plastic deformation parameters, and initial density of the powder, were measured. Multiple regression analysis demonstrated that the simulation parameters were significantly affected by process variables. The constructed DPC models were fed into the analysis using the finite element method (FEM), and the mechanical behavior of pharmaceutical powders during the tableting process was analyzed using the FEM. The results of this analysis revealed that the residual stress distribution of tablets increased with increasing X₄. Moreover, an interaction between X₂ and X₃ also had an effect on shear and the x-axial residual stress of tablets. Bayesian network analysis revealed causal relationships between the process variables, simulation parameters, residual stress distribution, and pharmaceutical responses of tablets. These results demonstrated the potential of the FEM as a tool to help improve our understanding of the residual stress of tablets and to optimize process variables, which not only affect tablet characteristics, but also are risks of causing tableting problems.

Various stresses, such as shear and axial stress, remain in tablets after the tableting process, because of the induction of elastic recovery. This type of stress, which is termed the residual stress distribution of the tablet, affects tablet characteristics and causes tableting problems. For instance, tablet failure, in particular capping, is more likely to be associated with an intensive shear band formed during the decompression stage.^1,2) Therefore, controlling the residual stress distribution of tablets is crucial in pharmaceutical design.

However, it is difficult to measure the residual stress distribution of tablets using analytical instruments. For this reason, the finite element method (FEM), in which the powder is modeled using the Drucker–Prager cap (DPC) model, is used to estimate the residual stress distribution of tablets.^3,4) The FEM is a numerical analytical method that is well established for the modeling of the deformation of powders in various industries, such as compaction in ceramic industries, and for the analysis of pharmaceutical-powder compaction.^5,6) In the FEM, powders are modeled as continuum media and the compaction behavior is analyzed by solving boundary value problems.

The DPC model is one of the continuum mechanical models, in which the powder is considered as a porous medium. The DPC model can represent the densification and hardening of the powder, as well as the interparticle friction.^7,8) This model can also reasonably represent both the shear failure and the plastic yielding of the powder, and can be readily characterized via experiments on real powders. Therefore, the DPC model is used frequently to analyze the strain, relative-density changes, and stress distribution of tablets during the tableting process. Finally, the DPC model is characterized by parameters that are mainly related to the internal friction angle, elastic modules, and plastic deformation, among others.

In general, the full calibration of the DPC model requires triaxial, hydrostatic compression, and proportional loading tests. These tests are used commonly for metallurgical powders.^9,10) Because pharmaceutical powders are very soft and loosely packed, the application of triaxial cells in pharmaceutical materials is difficult. To solve this problem, a novel method adapted for pharmaceutical powders was proposed.¹⁾ In this method, DPC model parameters are obtained by measuring axial upper/lower punch forces and displacements, as well as the radial die-wall pressure, during the tableting process.

Several studies have reported the calculation of the residual stress distribution of tablets using the FEM, in which the powder is modeled using the DPC model.^11–13) For instance, Kadiri and Michrafy reported that the stress distribution of tablets was affected by punch geometry.¹⁴⁾ It has been shown that the density distribution patterns are comparable with the experimental results of others.^15–18) Conversely, improvement of the DPC model was also studied.¹⁹⁾ In those studies, the DPC model parameters were considered as a function of the relative density of powders.

We have demonstrated that the residual stress distribution of tablets was affected by formulation and was closely related to the characteristics of the tablets, such as tensile strength and disintegration time.²⁰⁾ However, a few issues remained unclear. The first was that a placebo powder and one other type of powder (such as microcrystal cellulose or lactose) were used to model formulation in many studies,^1,5,19) even though the actual pharmaceutical products included one or more active pharmaceutical ingredients and many different excipients. Second, the effect of variables of the manufacturing process on DPC model parameters and residual stress distribution remains obscure, although the impact of the compression force has been well reported. Finally, the relationships between the DPC model parameters and pharmaceutical responses have not been elucidated. As the DPC model parameters express the physical characteristics of tablets, pharmaceutical characteristics, such as disintegration time, and dissolution property, may be predicted based on those parameters.

Therefore, the purpose of this paper was to investigate the problems described above. To achieve this purpose, we prepared 27 types of theophylline tablets according to the experimental design and measured the DPC model parameters of each tablet. Subsequently, we estimated the residual stress distribution of tablets using FEM and measured pharmaceutical characteristics of the tablets. Finally, we analyzed the data obtained using statistical methods such as multiple linear regression analysis (MRA) and Bayesian network (BN) analysis, to clarify the relationships between the process variables, DPC model parameters, residual stress distribution, and pharmaceutical characteristics of the tablets.

Experimental

Materials

Theophylline (JP grade) was purchased from Hachidai Pharmaceutical Co., Ltd. (Osaka, Japan). Lactose (LAC; Tablettose 80, Meggle Japan Co., Ltd., Tokyo, Japan), cornstarch (CS; Graflow M, Nippon Starch Chemical Co., Ltd., Osaka, Japan), and microcrystalline cellulose (MCC; Ceolus PH-101, Asahi Kasei Chemicals Co., Ltd., Tokyo, Japan) were purchased. Magnesium stearate (Mg-St) was purchased from Wako Pure Chemical Industries, Ltd. (Osaka, Japan). Polyvinylpyrrolidone K30 (PVP; Kollidon K30, BASF, Ludwigshafen, Germany) was a gift from the Nippon Shokubai Company, Ltd. (Osaka, Japan).

Experimental Design

To develop systematic model formulations, 27 types of tablets containing theophylline, LAC, CS, MCC, PVP, and Mg-St were prepared according to a Box–Behnken design (Table 1). Water amount (X₁), kneading time (X₂), lubricant-mixing time (X₃), and compression force (X₄) were selected as the process variables. Process conditions were fixed according to the flow chart presented in Table 2. All ingredients were dried at 75°C for 12 h. The ingredients were accurately weighed according to the experimental formulations, and all ingredients, with the exception of Mg-St, were blended using a mixer (KM4005; De’Longhi, Treviso, Italy) for 1 min. Distilled water was added as a granulation liquid and the mixture was kneaded (impeller speed, 470 rpm; processing time, 1–9 min). After the granulation process, the granules were sieved through a 5.8 mm mesh. The granules were dried at 75°C for 50 min and forcibly sieved through a 1.4 mm mesh. Mg-St was added to the granules and the mixture was blended using a V blender (S-3; Tsutsui Rikagaku Kikai, Tokyo, Japan) for 5–55 min at 50 rpm. The final blended product was compressed into flat-faced tablets (200 mg, 8 mm in diameter) using a tableting machine (AUTOTAB-500; Ichihashi-Seiki Company, Ltd., Kyoto, Japan). All the formulations were composed of 40 mg of theophylline, 86 mg of LAC, 37 mg of CS, 31 mg of MCC, 4 mg of PVP, and 2 mg of Mg-St.

Table 1. Box–Behnken Design of the Four Process Variables

Rp.	Water amount (%)	Kneading time (min)	Lubricant-mixing time (min)	Compression force (kN)
1	15	1	30	8
2	15	9	30	8
3	35	1	30	8
4	35	9	30	8
5	25	5	5	6
6	25	5	5	10
7	25	5	55	6
8	25	5	55	10
9	15	5	30	6
10	15	5	30	10
11	35	5	30	6
12	35	5	30	10
13	25	1	5	8
14	25	1	55	8
15	25	9	5	8
16	25	9	55	8
17	15	5	5	8
18	15	5	55	8
19	35	5	5	8
20	35	5	55	8
21	25	1	30	6
22	25	1	30	10
23	25	9	30	6
24	25	9	30	10
25	25	5	30	8
26	25	5	30	8
27	25	5	30	8

Table 2. Flow Chart of the Granulation Process, Simulation Parameters (α^y, E, ν, W₁^c, D₁^c, D₂^c, and D_i), Process Variables (X₁, X₂, X₃, and X₄), and Pharmaceutical Responses (TS, DT, and D₁₀)

Process	Condition	Process variables	Simulation parameters and pharmaceutical responses obtained experimentally
Drying	75°C, 12 h
Sieving	500 µm mesh
Blending	1 min, 470 rpm
Granulation	15–35%, 1–9 min, 470 rpm	Water amount (X1), kneading time (X2)
Sieving	5.4 mm mesh
Drying	75°C, 50 min
Sieving	1.4 mm mesh
Blending	V blender, 5–55 min, Mg-St 1%	Lubricant-mixing time (X3)	Internal friction angle (αy)
Tableting	6–10 kN, flat-faced tablet, 8 mm diameter	Compression force (X4)	Young’s modulus (E), Poisson rate (ν), plastic deformation parameters (W1c, D1c, D2c), initial density (Di)
Testing			Tensile strength (TS), disintegration time (DT), percent dissolved at 10 and 30 min (D10 and D30)

The simulation parameters were measured via direct shear test and analysis of compaction behavior.

Simulated Conditions of the Tableting Process Using FEM

Because the compaction of cylindrical tablets is carried out in an axisymmetric case, it can be analyzed using a two-dimensional FEM. Figure 1 shows the scheme of the FEM. The powder was considered as a DPC model. The die wall and upper punches were modeled as rigid bodies. The interaction between the powder, die wall, and upper punch was modeled, and the friction in the contacts was set to 0.1. The nodes on the symmetry axis were restricted to move only horizontally, and the nodes at the bottom boundaries were allowed to move only vertically. The upper punch could move vertically with compression.

Fig. 1. A Typical Finite Element Model for Modeling the Compaction of Flat-Faced Tablets

The powder was modeled using the DPC model. An axisymmetric two-dimensional model (right half) was used.

DPC Model

The DPC model is one of the yield surface models and can represent the plastic deformation, elastic deformation, and internal friction of powders. The constitutive equations of the model are not listed here, as they can be found in a previous paper.²⁰⁾ The DPC model parameters that have to be measured from the experiments are mainly Young’s modulus (E), Poisson rate (ν), internal friction angle (α^y), and three plastic deformation parameters (W₁^c, D₁^c, and D₂^c). The methods for measuring each parameter are described in the following sections.

Measurement of the Failure Envelope Using the Direct Shear Test

A direct shear tester (NS-V100; Nanoseeds Corporation, Gifu, Japan) was used to measure the failure envelope and estimate the internal friction angle (α^y). Before measurement, all ingredients were dried at 75°C for 12 h. The powder (2.5 g) was added into the shear cell. Shear stresses were then measured when the powder bed was compressed at 20, 40, and 60 N, respectively. The data observed were plotted in the axial shear (σ)–stress (τ) plane. The linear approximation method was used to estimate the failure envelope. The failure envelope was measured for three powders of each formulation. The details of the failure envelope are given in a previous paper.²⁰⁾

Measurement of Elastic Modules, Plastic Deformation Parameters, and Initial Density Using Analysis of Compaction Behavior

Elastic modules (Young’s modulus (E) and Poisson rate (ν)) and plastic deformation parameters (W₁^c, D₁^c, and D₂^c) were estimated based on references.^1,20) Powder compaction tests were conducted using an instrumented hydraulic press (TK-TB20KN; TOKUSHU KEISOKU Co., Ltd., Kanagawa, Japan). The axial upper/lower punch forces and displacements, and the radial die-wall pressure were measured during compaction. To obtain accurate measurements of the parameters, compression and decompression speeds were set at a relatively low speed (compression speed, 1 mm/s; decompression speed, 1 mm/s) and a large amount of sample powder was used (approximately 350 mg). The relative packing density at 2 MPa was measured based on the volume of the gap between the upper and lower punches. Moreover, the initial density (D_i) of the powder bed was calculated from this result. The E, ν, and D_i values were obtained from an average of three tablets, respectively.

A statistical method was applied to estimate the three DPC model parameters (W₁^c, D₁^c, and D₂^c) using volume change as an indicator. The W₁^c, D₁^c, and D₂^c were assigned to a Box–Behnken design, and 13 kinds of DPC models were constructed. The range of these parameters was determined using a preliminary test. The α^y value was obtained from the direct shear test measurements. R_c^y=0.6, R_t^y=10, X_i=10, σ_I=0.91, and φ=1 were set as arbitrary values that do not cause divergence problems in FEM analysis. The tableting process was simulated using the FEM and the thickness of tablets during compaction was measured at 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 kN. The data observed were modeled by RSM-S, and DPC model parameters that were correlated with the experiment were estimated.

Evaluation of the Pharmaceutical Responses of Tablets

The hardness of the tablets was determined using a tablet-hardness tester (Portable checker PC-30; Okada Seiko, Tokyo, Japan). TS was calculated as:

  

(1)

where F is the maximal diametrical crushing force and d and t are the diameter and thickness of the tablet, respectively. TS values were measured for three tablets of each formulation.

The disintegration test was performed according to the JP16 disintegration test for tablets using a disintegration tester (NT-20H; Toyama Sangyo Co., Ltd., Osaka, Japan) and water (as a test medium) at 37°C. The DT was defined as the interval that was required for the complete disappearance of a tablet or its particles from the tester net. The DT was measured for three tablets of each formulation.

The dissolution test was performed according to the JP16 dissolution test No. 2 (the paddle method) at 50 rpm (NTR-6100A; Toyama Sangyo Co., Ltd.). The dissolution medium used was 900 mL of distilled water at 37±0.5°C. The samples were collected and filtered after 1, 4, 10, 20, 30, 45, and 60 min. The concentration of theophylline was measured spectrophotometrically at 271 nm using a Jasco Ubest-30 spectrophotometer (Japan Spectroscopic Co., Ltd., Tokyo, Japan). The dissolution rates of three tablets of each formulation were measured.

Evaluation of the Homogeneity of Residual Stress Distribution

The standard deviation (S.D.) was used as an indicator of homogeneity. S.D. is defined as:

  

(2)

where x̅, x_i, and n indicate the average of the stress distribution of each formulation, each stress point obtained by FEM, and the number of stress points obtained by FEM, respectively. Because n was determined by the number of meshes used in the FEM, n=333 was used in this study.

Modeling of the Causal Relationships between Factors

To evaluate the causal relationships between factors, we applied a nonlinear response surface method that incorporated a thin plate spline interpolation (RSM-S), which has been used to determine acceptable formulations of pharmaceutical compounds. Using RSM-S, we can easily understand the nonlinear relationships between causal factors and response variables and estimate a robust optimal solution.²¹⁾

MRA was performed to clarify the impact of process variables on the simulation parameters or residual stress distribution of tablets. A forward selection method based on stepwise selection was applied to the selection of causal factors, and factors with p values >0.25 were successively excluded from the analysis. In general, p>0.25 is used as a judging standard for excluding factors. Lowering the criterion is increasing the risk of involving meaningless factors, meanwhile significant factors are possibly eliminated with tightening the criterion.

Bayesian networks (BN) were used to construct the probabilistic graphical model among the variables of the manufacturing process, DPC model parameters, residual stress distribution of tablets, and tablet properties, and to estimate conditional independencies. BN are popular in statistics, machine learning, and artificial intelligence, and are mathematically strict and intuitively understandable. They represent the probabilistic relationships between random variables by using a directed acyclic graph and a set of conditional probability distributions.²²⁾ The construction of BNs requires candidates for the parent node (explanatory variable). In this study, a three-layered BN model was constructed, i.e., the process variables were assigned as the parent nodes of the DPC model parameters, and the DPC model parameters were assigned as the parent nodes of the stress distribution of tablets and of the pharmaceutical responses of tablets. The factors were discretized into three categories (low, medium, and high) using the K-means method.

Computer Programs

The FEM analysis of the tableting process was performed using ANSYS® 14.5 (ANSYS Inc., Canonsburg, PA, U.S.A.). The RSM-S was performed using dataNESIA® version 3.2 (Azbil Corporation, Tokyo, Japan). The MRA was performed using JMP version 8 (SAS Institute Inc., Cary, NC, U.S.A.). The BN model was constructed using BayoNet System software, version 6.0 (Mathematical Systems Inc., Tokyo, Japan).

Results and Discussion

Evaluation of the Pharmaceutical Responses of Tablets

Measured data, such as tensile strength (TS), disintegration time (DT), and dissolution property, are summarized in Table 3. Percent dissolved at 10 min (D₁₀) was selected as a typical dissolution property of tablets at the early stage, because the results of a dissolution test revealed the presence of a great difference in each sampling point at 10 min. Although percent dissolved at 30 min (D₃₀) was selected as a late-stage dissolution property, D₃₀ was only used in BN analysis.

Table 3. Pharmaceutical Responses of All Formulations

Rp.	Tensile strength (MPa)	Disintegration time (s)	Percent dissolved at 10 min	Percent dissolved at 30 min
1	1.45±0.12	83±8	78.4±5.0	94.3±2.6
2	1.66±0.11	123±27	64.9±1.4	92.4±1.8
3	1.78±0.06	1730±44	37.4±6.0	70.4±4.2
4	0.91±0.09	418±106	53.8±16.2	91.7±2.0
5	0.94±0.00	60±12	74.7±2.3	95.2±0.5
6	1.83±0.01	920±45	86.7±5.8	99.7±0.4
7	0.78±0.09	66±3	46.9±5.0	86.6±5.6
8	1.59±0.07	434±29	38.4±2.5	86.9±6.7
9	0.86±0.04	149±29	74.6±3.1	97.8±0.7
10	1.62±0.15	186±63	82.3±6.8	96.0±0.2
11	0.68±0.09	205±23	49.1±7.2	86.9±3.3
12	1.35±0.03	1620±107	16.1±0.7	41.3±3.0
13	1.89±0.09	110±22	87.0±0.7	104.6±0.7
14	1.50±0.08	551±117	79.2±14.8	103.7±4.5
15	1.26±0.18	56±8	75.1±6.0	91.9±1.9
16	0.95±0.11	80±20	62.9±6.3	92.1±1.8
17	1.64±0.06	71±19	64.4±2.8	95.5±3.3
18	1.30±0.26	85±17	71.3±3.9	95.5±1.6
19	1.22±0.07	512±68	28.3±2.1	57.0±2.5
20	0.99±0.11	1125±78	22.2±2.5	48.6±4.8
21	1.12±0.07	56±6	76.4±3.0	92.5±1.3
22	1.83±0.19	557±76	50.4±4.4	94.1±1.3
23	0.66±0.08	103±11	69.8±2.3	90.2±1.3
24	1.36±0.18	462±20	73.5±1.5	90.0±0.3
25	1.25±0.17	344±69	60.3±10.7	90.1±5.4
26	1.21±0.07	284±88	72.7±1.5	95.1±1.8
27	1.21±0.09	249±47	53.8±5.9	91.6±1.6

Each data point is the mean±S.D. (n=3).

The response surfaces for TS, DT, and D₁₀ were estimated using RSM-S based on the original data set. The accuracy of the response surfaces was evaluated via leave-one-out cross-validation (LOOCV), which revealed that the correlation coefficients for TS, DT, and D₁₀ were sufficiently high (0.934, 0.959, and 0.888, respectively). Multiple regression analysis was also performed. The coefficient of determination, which was adjusted using degrees of freedom (R**²) and is an indicator of the fit of each linear regression equation, was estimated. The R**² values for TS, DT, and D₁₀ were 0.827, 0.765, and 0.714, respectively, resulting in poor estimations.

As shown in Fig. 2, the results regarding response surfaces revealed that the TS increased as the kneading time (X₂) and lubricant-mixing time (X₃) decreased and the compression force (X₄) increased. The DT increased with increasing water amount (X₁), X₃, and X₄, and decreasing X₂. The D₁₀ was strongly affected by X₁ and X₃, and increased as X₁ and X₃ decreased.

Fig. 2. Response Surfaces for Tensile Strength (a), Disintegration Time (b), and Percent Dissolved at 10 min (c) as a Function of Water Amount, Kneading Time, Blending Time, and Compression Force

Eighty-one data points obtained by experimental design were analyzed based on RSM-S. LOOCV results showed high prediction ability in all models (more than at least 0.88).

Granule properties, such as mean particle size, porosity, and specific surface area, are responsible for tablet characteristics. For instance, Ohno et al. reported that both the 50% pore diameter and the 50% particle diameter have a strong influence on the dissolution property.²³⁾ Moreover, the 50% pore diameter decreases with increasing water amount. Consequently, a greater X₁ decreased DT and D₁₀ because excessive water induced a decrease of 50% in pore diameter in this study. Mg-St reduces the interparticle binding strength, DT, and dissolution property, because Mg-St has hydrophobic and flatting characteristics. Therefore, excessive mixing induces a delaying effect regarding the reduction of the hardness of the tablet and the release rate of the drug.

Effect of Process Variables on Simulation Parameters

The results of the measurement of Young’s modulus (E), Poisson rate (ν), internal friction angle (α^y), plastic deformation (W₁^c, D₁^c, and D₂^c), and initial density (D_i) of each formulation are summarized in Table 4. To understand visually the effect of process variables on simulation parameters, response surfaces for E and ν were estimated using RSM-S, as shown in Fig. 3. The accuracy of the response surfaces was evaluated via LOOCV, which revealed that the correlation coefficients for E and ν were sufficiently high (0.710 and 0.729, respectively). To evaluate the contribution of each factor to simulation parameters, multiple regression analysis (MRA) was also performed. The results are summarized in Table 5.

Table 4. Simulation Parameters of All Formulations

Rp.	Elastic modules		Friction	Bulkiness	Plastic deformation
	Young’s modules (GPa)	Poisson rate	Internal friction angle (°)	Initial density (mg/mm3)	W1c	D1c	D2c
1	7.76±0.41	0.1255±0.0029	29.1±4.9	0.778±0.014	0.511	0.0181	5.72E–09
2	8.84±0.63	0.1259±0.0031	22.6±1.6	0.745±0.012	0.564	0.0172	6.52E–09
3	8.68±0.55	0.1221±0.0031	24.5±2.6	0.722±0.027	0.593	0.017	5.29E–09
4	8.19±0.45	0.1205±0.0017	23.3±2.0	0.760±0.029	0.549	0.0153	1.70E–09
5	7.40±0.88	0.1178±0.0034	22.0±3.4	0.741±0.008	0.545	0.0178	2.57E–09
6	8.79±0.87	0.1257±0.0045	22.0±3.4	0.737±0.021	0.557	0.0172	5.17E–09
7	7.40±0.83	0.1164±0.0022	19.4±1.6	0.749±0.019	0.512	0.0183	4.08E–09
8	8.37±0.52	0.1260±0.0049	19.4±1.6	0.754±0.010	0.538	0.0174	4.83E–09
9	6.51±0.25	0.1198±0.0015	23.2±5.2	0.742±0.019	0.53	0.0182	8.18E–09
10	8.12±0.05	0.1294±0.0032	23.2±5.2	0.758±0.015	0.535	0.0171	5.03E–09
11	6.93±0.72	0.1170±0.0012	22.1±3.5	0.740±0.030	0.541	0.0184	2.50E–09
12	8.89±0.04	0.1282±0.0032	22.1±3.5	0.727±0.006	0.6	0.0172	5.72E–09
13	8.66±0.06	0.1198±0.0011	20.7±3.9	0.705±0.012	0.616	0.017	5.50E–09
14	7.49±0.41	0.1223±0.0019	20.4±4.7	0.739±0.004	0.565	0.0171	5.76E–09
15	7.76±0.42	0.1227±0.0012	16.8±1.6	0.727±0.020	0.597	0.0155	5.18E–09
16	8.72±0.54	0.1197±0.0036	18.8±1.7	0.703±0.018	0.617	0.0168	5.31E–09
17	8.69±0.02	0.1253±0.0012	22.1±2.6	0.721±0.019	0.601	0.0161	5.47E–09
18	7.64±0.62	0.1234±0.0005	23.2±3.0	0.754±0.017	0.537	0.0176	6.91E–09
19	7.92±0.49	0.1241±0.0024	24.4±2.0	0.739±0.010	0.564	0.0171	6.00E–09
20	8.38±0.49	0.1212±0.0022	23.1±2.3	0.694±0.014	0.626	0.0177	5.50E–09
21	7.40±0.53	0.1165±0.0028	19.0±2.0	0.726±0.023	0.553	0.0185	4.23E–09
22	9.84±0.06	0.1239±0.0030	19.0±2.0	0.725±0.011	0.6	0.0172	6.02E–09
23	7.43±0.79	0.1128±0.0027	19.3±2.2	0.741±0.023	0.543	0.015	5.50E–09
24	9.19±0.52	0.1245±0.0008	19.3±2.2	0.706±0.036	0.604	0.0201	8.41E–09
25	7.90±0.67	0.1201±0.0020	18.8±1.6	0.737±0.024	0.57	0.0171	6.17E–09
26	8.05±0.22	0.1181±0.0005	22.3±1.9	0.716±0.015	0.609	0.017	5.65E–09
27	8.50±0.29	0.1209±0.0006	18.6±1.2	0.702±0.035	0.612	0.0166	4.83E–09

Each data point (Young’s modulus, Poisson rate, internal friction angle, and initial density) is the mean±S.D. (n=3).

Fig. 3. Response Surfaces for Simulation Parameters: Young’s Modulus (a and b) and Poisson Rate (c and d)

Eighty-one data points obtained by experimental design were modeled as a function of water amount, kneading time, lubricant-mixing time, and compression force using RSM-S. Background factors were set to intermediate values in all cases.

Table 5. Prediction Accuracy and p Value of Multiple Regression Analysis

	Elastic modules		Internal friction	Bulkiness	Plastic deformation
	E	ν	αy	Di	W1c	D1c	D2c
R**2	0.777	0.922	0.648	0.388	0.153	0.052	0.151
RMSE	0.349	0.001	1.517	0.016	0.032	0.001	0.000
Water amount (X1)	0.253	0.0005*	0.482	0.0479*	0.089	—	0.0257*
Kneading time (X2)	0.819	0.300	0.0249*	0.815	—	0.170	—
Lubricant-mixing time (X3)	0.331	0.108	—	0.680	0.447	0.246	—
Compression force (X4)	<0.0001*	<0.0001*	—	—	—	—	—
X12	—	<0.0001*	<0.0001*	0.051	—	—	—
X22	0.0149*	—	—	—	—	—	—
X32	—	0.0357*	—	—	—	—	—
X42	—	—	—	—	—	—	—
X1×X2	0.0378*	0.0222*	—	0.0371*	—	—	—
X1×X3	0.0440*	—	0.094	0.0235*	0.059	—	—
X1×X4	—	—	—	—	—	—	—
X2×X3	0.0067*	—	—	0.083	—	—	—
X2×X4	—	0.069	—	—	—	—	—
X3×X4	—	—	—	—	—	—	—

Each simulation parameter was predicted based on the four process variables. A forward selection method based on stepwise selection was applied to the selection of causal factors, and factors with p values >0.25 were successively excluded from the analysis.

E was significantly affected by X₄, X₂², X₁×X₂, X₁×X₃, and X₂×X₃, as shown in Table 4. In particular, the effects of X₄ and X₂×X₃ were strong. The analysis of response surfaces showed that a higher X₄ led to an increased E. E was slightly decreased with decreasing X₃ (Fig. 3b). Although E increased as X₂ decreased when X₁ was high, E increased as X₂ increased when X₁ was low (Fig. 3a). A higher E indicates tablets that are more resistant to elastic deformation, because the E of materials represents the stress that is required to produce a corresponding unit of strain. The results of the MRA showed the ν was driven by X₁, X₄, X₁², X₃², and X₁×X₂. The effect of X₁ on ν was comparable to that of X₄. Response surfaces showed that ν increased with decreasing X₁ and increasing X₄ (Figs. 3c, d). The impact of compression force on E and ν observed here was consistent with the findings of a previous report.¹⁾

The α^y was affected by X₁, X₂, X₁², and X₁×X₃ (Table 5), with the impact of X₂ and X₁×X₃ being particularly strong. Rp.1 showed the highest value among the 27 kinds of test formulations (Table 4). Because Rp. 1 was set to have a small X₁ and short X₂, granulation might not have progressed much. This difference might affect the experimental results. D_i was affected by X₁, X₂, X₃, X₁², X₁×X₂, X₁×X₃, and X₂×X₃. More importantly, X₁, X₁×X₃, and X₂×X₃ contributed strongly to changes in D_i. In general, a higher granule size yields greater bulk, and granule size is affected by X₁ and X₂. Conversely, as X₃ affects interparticle friction and agglomeration, this effect might contribute to D_i. The model of W₁^c, D₁^c, and D₂^c based on MRA yielded poor estimation value and the effect of factors could not be evaluated well. A possible explanation for this low prediction ability is strong nonlinearity, because plastic deformation of the DPC model was correlated with experimental values and could be modeled by a BN, as shown in the following section.

Comparison of the FEM with Experimental Results

The material properties that were obtained from experiments were fed into the ANSYS software, in which the DPC model was implemented. The mean values obtained in experiments were used as the DPC model parameters, as shown in Table 4. In this study, appropriate values of the parameters R_c^y, R_t^y, X_i, σ_I, and φ were sought in an arbitrary manner when they could not be measured based on the data, because these parameters have little effect on plastic deformation and stress distribution.²⁰⁾

In many cases, the volume change of the powder during the tableting process was used to compare FEM with experimental results, because there is no concise procedure that enables the simple visualization of the residual stress distribution of tablets.^1,11) Figure 4a shows a typical example of the variation of axial stress with axial strain during compaction. This result showed that the FEM and experimental curves were very close. To evaluate all formulations, the axial strain of the compact obtained by FEM was plotted as a function of the axial strain of tablets obtained based on experimental results, and the linear approximation method was applied. Subsequently, correlation coefficient, gradient, and intercept of liner were calculated. The results of tablets and compacts for each compression force are shown in Fig. 4b and Table 6, respectively. This result showed high correlation coefficients in many cases, although a slightly poor estimation was observed at a low compression force (1 and 2 kN). Moreover, the values of gradient and intercept were relatively close to one and zero, respectively, indicating good correlations. These results suggest that FEM results corresponded well to experimental findings.

Fig. 4. Comparison of FEM and Experimental Results: (a) Variation of Axial Stress with Axial Strain during Compaction; (b) Axial Strain of Each Tablet (Rp. 1–27)

In all cases, the experimental values obtained were in good agreement with the FEM results.

Table 6. Correlation Coefficient, Gradient, and Intercept of a Straight Line That Represents the Relationship between the Axial Strain of FEM and Experimental Results at Each Compression Force

Compression force (kN)	1	2	3	4	5	6	7	8	9	10
Correlation coefficient	0.872	0.936	0.926	0.935	0.931	0.953	0.971	0.974	0.961	0.941
Gradient of the line	0.577	0.737	0.956	1.073	1.138	1.141	1.045	1.020	0.942	1.001
Intercept of the line	0.086	0.099	0.035	−0.015	−0.051	−0.061	−0.023	−0.016	0.014	−0.021

Impact of Process Variables on Residual Stress Distribution of Tablets

The typical residual stress distributions of tablets that showed the most homogeneous and inhomogeneous formulations are shown in Fig. 5. We evaluated the effect only of X₁, X₂, and X₃, because it is well known that the residual stress increases with increasing X₄. Consequently, tablets compressed at 8 kN were used for the comparison of residual stress distribution.

Fig. 5. Two-Dimensional Maps for the Residual Shear and x-Axial and y-Axial Stress Distribution of Tablets Estimated Using the FEM

As typical examples, the most homogeneous and inhomogeneous formulations are shown here. The process variables of each formulation are shown in Table 2.

The results of the FEM analysis revealed that the trend of the residual stress distribution was similar among model formulations. An intensive shear region was recognized from the top edge to the mid center. In this shear region, the shear stress changed from positive to negative, indicating a change in the direction of the shear stress. Regarding x-axial stress (σ_x), the top corner exhibited the highest value, whereas the bottom corner exhibited the lowest value. It was clear that the absolute value of y-axial stress (σ_y) at the top was generally lower than that observed at the bottom. Moreover, the σ_y value at the top had the stress extending upward, whereas the one at the bottom had the stress extending downward. The overall trend was similar to that reported previously.²¹⁾

Although the residual stress distribution was similar among model formulations, the ranges of stress values were significantly different. The MRA revealed that the residual stress distribution of tablets was significantly affected by X₄, as it increased with increasing X₄. Moreover, an interaction between X₂ and X₃ also affected shear and the x-axial residual stress of tablets significantly. To evaluate the impact of process variables on residual stress distribution, response surfaces for S.D. of each stress distribution as a function of four process variables were estimated using RSM-S (Fig. 6). A high X₁ and short X₂ resulted in low S.D. values, indicating a uniform stress distribution (Figs. 6a, c, e). Conversely, higher X₃ and X₄ yielded a more inhomogeneous stress distribution in many cases (Figs. 6b, d, f). There was an optimum blending time of the lubricant, because overblending led to agglomeration and a decrease of interparticle bonding. This phenomenon might cause higher residual stress.

Fig. 6. Response Surfaces for Simulation Parameters: S.D. of Shear Stress (a and b), S.D. of σ_x Stress (c and d), and S.D. of σ_y Stress (e and f), as a Function of Water Amount, Kneading Time, Blending Time, and Compression Force

Twenty-seven data points obtained by FEM were analyzed based on RSM-S. LOOCV results showed high prediction ability in all models. Background factors were set to intermediate values in all cases.

We performed a successful quantitative prediction of tablet characteristics, as well as residual stress distribution, using FEM and RSM-S. Consequently, we were able to seek optimum process variables that were considered not only as tablet characteristics, but also as a risk of having tableting problems.

Construction of a BN Model

BNs efficiently implement the probabilistic inference algorithm, which estimates the probability distribution of arbitrary random variables in a model.^24,25) To analyze the latent structure among the process variables, FEM parameters, residual stress distribution, and pharmaceutical responses of tablets, BN models were constructed using Akaike’s information criteria (AIC), the K2 algorithm, and minimum description length (MDL) as the judging criteria. The distinctive probabilistic model was estimated by the edges, which represent conditional dependencies, and between the nodes, which represent the variables.

The internal structures of BN models differed slightly depending on the judging criteria used. Therefore, the optimal BN model was estimated using the indices of accuracy rate, precision, recall, and F-measure. In general, there is a trade-off between precision and recall, as greater precision decreases recall and vice versa. The F-measure is the harmonic mean of precision and recall and takes both measures into consideration. When simulation parameters, the S.D. of each stress, and the pharmaceutical responses of tablets were predicted based on process variables, all measures estimated using the model based on the K2 algorithm were greater than 0.77, indicating that the prediction ability of the probabilistic model is sufficiently high. In contrast, all measures estimated using the model based on AIC and MDL were lower than 0.66 (Table 7). These results led us to conclude that the BN model based on the K2 algorithm has an optimal structure and represents a significant latent structure.

Fig. 7. Bayesian Network Model of the Latent Structure among Process Parameters, Simulation Parameters, S.D. of Residual Stress Distribution, and Responses of Tablets Estimated Using the K2 Algorithm

Table 7. Four Standard Measures of Bayesian Network Models Based on Each Measurement Criterion

	Accuracy rate	Precision	Recall	F-Measure
K2 algorithm	0.829	0.800	0.777	0.772
Akaike’s information criteria (AIC)	0.656	0.593	0.567	0.549
Minimum description length (MDL)	0.598	0.386	0.494	0.421

The selection of simulation parameters as input values allowed the calculation of the accuracy rate of the other parameters. This result showed that the accuracy rate of pharmaceutical responses ranged from 0.852 to 0.864, indicating a high prediction ability. This finding suggests that the simulation parameters are closely related to pharmaceutical responses, such as TS, DT, D₁₀, and D₃₀. However, it might be difficult to predict pharmaceutical responses quantitatively based on simulation factors, because prediction ability decreased drastically when the factors were discretized into five categories. In addition, acceptable results were not obtained based on MRA.

Evaluation of the Relationships between Process Variables, Simulation Parameters, Stress Distribution, and Pharmaceutical Responses Based on the BN Model

The BN model constructed here is shown in Fig. 7. The use of a BN model allowed the thorough understanding of the relationships between the variables of the manufacturing process, DPC model parameters, residual stress distribution, and pharmaceutical responses of tablets. The BN model obtained clarified the causal relationships among factors. For instance, Young’s modulus (E) was affected by all process variables, and affected TS, DT, and S.D. of τ. The water amount (X₁) and compression force (X₄) had an impact on the Poisson rate (ν), and ν was correlated with TS, DT, and S.D. of τ and σ_y. X₁, kneading time (X₂), and blending time (X₃) had an effect on the α^y. Finally, α^y was associated with TS, DT, D₁₀, D₃₀, and S.D. of σ_x. These results were similar to the MRA results (Table 5).

To evaluate semi-quantitative dependence, the posterior probability of the causal factors was predicted based on the BN model constructed here. Four typical examples of the results obtained are shown in Fig. 8. First, the effect of elastic modules was evaluated. The X₄ and S.D. of τ were relatively low at a high α^y, small E, and small ν. Conversely, a high α^y, high E, and high ν yielded a low X₁ and a high X₄ and S.D. of τ (Fig. 8a). This result indicates that E and ν increased with decreasing X₁, increasing X₄, and increasing S.D. of τ. Therefore, we focused on the variation of α^y. The BN model inferred factors with low α^y, intermediate E, and intermediate ν for intermediate X₁, TS, and D₃₀ tablet values. Conversely, there was no remarkable difference in CPDs at high α^y, intermediate E, and intermediate ν (Fig. 8b). This result indicates that α^y is a governing factor in this condition. Third, we focused on the variation of D_i. The BN model inferred factors with low D_i, intermediate D₁^c, and intermediate D₂^c for tablet intermediate X₃, intermediate S.D. of τ, and intermediate S.D. of σ^x. However, there was no remarkable difference in CPDs at high D_i, intermediate D₁^c, and intermediate D₂^c (Fig. 8c). Finally, we estimated the critical factors of D₁₀. When tablets had a small α^y, high D_i, small W_c1, intermediate D₁^c, and intermediate D₂^c, the BN model inferred parameters with intermediate X₂, high X₃, low D₁₀, and high S.D. of σ_x; in contrast, when tablets had small α^y, intermediate D_i, intermediate W_c¹, high D₁^c, and intermediate D₂^c, the BN model inferred parameters with low X₂, intermediate X₃, high D₁₀, and low S.D. of σ_x (Fig. 8d). Although a higher W_c¹ resulted in an increase of volume change at all compression forces, a higher D_c¹ indicated a higher volume change at low compression forces. These results suggest that D₁₀ is closely related to plastic deformation and bulkiness, because the D₁₀ value was dramatically changed when D_i, W_c¹, and D_c¹ were selected as input values.

Fig. 8. Typical Example of the Conditional Probability Distributions (CPDs) of Causal Factors Inferred Using a Bayesian Network Model

(a) CPDs of X₁, X₄, and S.D. of τ predicted from α^y, E, and ν; (b) CPDs of X₁, TS, D₃₀, and S.D. of σ_x predicted from α^y, E, and ν; (e) CPDs of X₁ and S.D. of τ, σ_x, and σ_y predicted from D_i, D₁^c, and D₂^c; (d) CPDs of X₁, X₃, and S.D. of τ and σ_y predicted from α^y, D_i, W₁^c, D₁^c, and D₂^c.

Conclusion

In this study, we investigated the relationships between process variables, residual stress distribution of tablets, and pharmaceutical responses. A direct shear test and instrumented die were used to model the mechanical behavior of pharmaceutical powders. The FEM was used to simulate the compaction behavior of the powders and estimate the residual stress distribution of each formulation. Simulation parameters, such as E, ν, α^y, and D_i, were expressed as a function of the four process variables using MRA, and the impact of process variables on simulation parameters was evaluated. The FEM revealed that the residual stress distribution of tablets was affected by the four process variables. The BN model clarified the causal relationships between process variables, simulation parameters, pharmaceutical responses, and residual stress distribution of tablets. These results demonstrated that FEM is a useful tool to help improve our understanding of residual stress and optimize process variables that are considered not only as tablet characteristics, but also as a risk of having tableting problems.

Acknowledgment

This study was supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The authors thank Mr. Yuuki Ishida, Department of Pharmaceutics, Hoshi University, for performing the dissolution test.

References

• 1) Wu C. Y., Ruddy O. M., Bentham A. C., Hancock B. C., Best S. M., Elliott J. A., Powder Technol., 152, 107–117 (2005).
• 2) Wu C. Y., Hancock B. C., Mills A., Bentham A. C., Best S. M., Elliott J. A., Powder Technol., 181, 121–129 (2008).
• 3) Coube O., Riedel H., Powder Metallurgy, 43, 123–131 (2000).
• 4) Sinka I. C., Cunningham J. C., Zavaliangos A., Powder Technol., 133, 33–43 (2003).
• 5) Cunningham J. C., Sinka I. C., Zavaliangos A. J., J. Pharm. Sci., 93, 2022–2039 (2004).
• 6) Michrafy A., Diarra H., Dodds J. A., Michrafy M., Penazzi L., Powder Technol., 208, 417–422 (2011).
• 7) Aydin I., Briscoe B. J., Sanliturk K. Y., Powder Technol., 89, 239–254 (1996).
• 8) Aydin I., Briscoe B. J., Ozkan N., MRS Bull., 22, 45–51 (1997).
• 9) DiMaggio F. L., Sandler I., ASCE, 97, 935–950 (1971).
• 10) Häggblad H. A., Powder Technol., 67, 127–137 (1991).
• 11) Diarra H., Mazel V., Boillon A., Rehault L., Busignies V., Bureau S., Tchoreloff P., Powder Technol., 224, 233–240 (2012).
• 12) Sinka I. C., Cunningham J. C., Zavaliangos A., J. Pharm. Sci., 93, 2040–2053 (2004).
• 13) Michrafy A., Ringenbacher D., Tchoreloff P., Powder Technol., 127, 257–266 (2002).
• 14) Kadiri M. S., Michrafy A., Powder Technol., 239, 467–477 (2013).
• 15) Eiliazadeh B., Briscoe B. J., Sheng Y., Pitt K., Particul. Sci. Technol., 21, 303–316 (2003).
• 16) Laity P. R., Cameron R. E., Powder Technol., 192, 287–297 (2009).
• 17) Djemai A., Sinka I. C., Int. J. Pharm., 319, 55–62 (2006).
• 18) Sinka I. C., KONA Powder and Practice Journal, 25, 4–22 (2007).
• 19) Diarra H., Mazel V., Busignies V., Tchoreloff P., Int. J. Pharm., 453, 389–394 (2013).
• 20) Hayashi Y., Miura T., Shimada T., Onuki Y., Obata Y., Takayama K., J. Pharm. Sci., 102, 3678–3686 (2013).
• 21) Hayashi Y., Oshima E., Maeda J., Onuki Y., Obata Y., Takayama K., Chem. Pharm. Bull., 60, 1419–1425 (2012).
• 22) Cooper E. G., Herskovitz E., Mach. Learn., 9, 309–347 (1992).
• 23) Ohno I., Hasegawa S., Yada S., Kusai A., Moribe K., Yamamoto K., Int. J. Pharm., 338, 79–86 (2007).
• 24) Motomura Y., Synthesiology, 2, 1–11 (2009).
• 25) Bothner U., Milne S. E., Kenny G. N. C., Georgieff M., Schraag S., J. Clin. Monit. Comput., 17, 31–36 (2002).