Predicting Drug Release from 3D Printed Oral Medicines Based on the Surface Area to Volume Ratio of Tablet Geometry

3D printing offers the advantage of being able to modify dosage form geometry, which can be exploited to modify release characteristics. In this study, we investigated the influence of the surface area to volume ratio (SA/V) to change and predict release profiles of 3D printed dosage forms. Geometries with varying SA/V and dosages were designed and printed, and drug dissolution was investigated. Three drug substances were used: pramipexole, levodopa (both BCS I) and praziquantel (BCS II). Two polymers were chosen as matrix formers: polyvinyl alcohol (water-soluble) and ethylene vinyl acetate (inert). Drug release was characterized using the mean dissolution time (MDT) and established equations that describe complete dissolution curves were applied. Predictions were validated with previously un-printed dosage forms. Based on an identified MDT-SA/V correlation, the MDT can be predicted with a deviation of ≤5 min for a given SA/V. Using correlations of fit parameters and SA/V, RMSEP values of 0.6–2.8% and 1.6–3.4% were obtained for the BCS I formulations and RMSEP values of 1.0–3.8% were obtained for the BCS II formulation, indicating accurate prediction over a wide range of dissolution profiles. With this approach, MDT and release profiles of dosage forms with a given SA/V can be precisely predicted without performing dissolution tests and vice versa, the required SA/V can be predicted for a desired release profile.


Introduction
Personalized pharmaceutical therapies are increasingly moving into focus to match the individual needs of patients [1][2][3]. There is a constantly improving understanding of the organism and metabolism of children and elderly patients, and it is becoming increasingly clear that drug substances may be metabolized differently compared to adults, especially if comorbidities are present. The efficacy of the active pharmaceutical ingredient (API) varies depending on age, body type, gender, health, and the individual metabolic state [4][5][6]. To adjust the dose of medicine for patients based on these dependencies, several approaches are available. Liquid or semi-solid preparations can be measured and dosed individually with the help of dosing aids, such as spoons, syringes, or cups. With solid dosage forms, it becomes more difficult to provide exact doses. Some tablets can be split into halves or quarters by hand or with the help of a tablet cutter, but these methods lead to inaccurate dosing [7,8]. Tablets with modified release especially demonstrate a major challenge if the release-determining factor is controlled by a functional coating, which loses its integrity when a tablet is split or mortared [9]. The requirements for dose variability and the constraining demand for an equivalent release profile despite the change in dose strength cannot be addressed industrially. The small batches needed for personalized medicine are not technically and economically feasible to manufacture on large-scale equipment. To provide patients with customized medicines, 3D printing to manufacture individualized medicine has been investigated [10][11][12][13][14]. Most commercial 3D printers are not only inexpensive, but also easy to operate and able to produce small batches for patients on demand. With individually manufactured dosage forms, it is possible to adapt the dosage of the API and the release kinetics of the API to the patient. In this way, a wide range of therapeutic regimens can be covered and side effects reduced, as the most effective dose with the least adverse effects can be administered precisely [15]. In the future, the production of these individual batches could be realized by community pharmacies, hospitals and regional manufacturing hubs.
In fused deposition modeling (FDM) 3D printing, intermediates containing the API, so-called filaments, are produced by hot-melt extrusion [3,[16][17][18]. Filaments have to meet certain mechanical properties in order to enable the subsequent printing process [16,19,20]. The drug-loaded filament is fed into the print head, molten in the nozzle, and deposited on a heated print bed. Because of these two heat-intensive steps, FDM 3D printing is limited in terms of the APIs that can be used. The APIs should be thermostable and only decompose at high temperatures. Although it is possible to work with polymers that melt at low temperatures, this severely limits the choice of polymers [21]. The desired object is designed in advance using computer-aided design (CAD) software and converted via a slicer software into a machine language, the G-code. This G-code determines the temperature, movement, and the speed of the print head and print bed. The design of the object should consider the volume, which determines the dosage depending on the drugload of the filament. There are two opportunities to influence the drug release properties of 3D printed objects: changing the composition of the filaments and adapting the surface area to volume (SA/V) ratio [22][23][24][25]. Since printing technologies allow one to manufacture a wide range of geometries without further development efforts, changing the SA/V ratio is the preferred option. It has been shown that the API release is faster, when the SA/V ratio is higher and slower when the ratio is smaller [22,23,26]. The absolute volume and the selected geometry were found to be negligible [23]. In order to reduce the elaborate production of filaments, there is a need to achieve different release profiles with only one filament.
To enable an individual dosage form with a defined dissolution profile, one must understand how the underlying formulation performs and according to which kinetic the API is released from the polymer matrix. Release profiles can be fitted using mathematical equations that consider the dissolution behaviour of the API and the polymer, as well as the diffusion pathways the API and the medium must overcome [27][28][29]. So far, various approaches have been used to try to predict the drug release from solid dosage forms [30][31][32][33][34][35][36][37][38][39][40][41]. Korte et al. investigated the approach of changing the dose of the 3D geometry via the percent of infill and predicted the release kinetics depending on the infill density [42]. Likewise, in another study, an infill in the form of honeycomb was modified, thereby also changing the drug release profile and predicting the resulting release profiles [43]. Artificial neural networks (ANN) are also being implemented to identify the influence of formulation and process parameters on the release behaviour and to improve predictions [28,34,[44][45][46][47]. For example, Novák investigated the influence of varying infills and resulting tablet porosities on drug release using an ANN. This resulted in an in-silico design method of infill variations for 3D printed dosage geometries [48].
In this work, two different approaches were applied to predict dissolution characteristics and complete dissolution profiles of 3D printed geometries. This should ensure that individual batches can be produced in hospitals and community pharmacies without having to spend a lot of material, time and money on release testing. Firstly, the release behaviour of various geometries should be predicted via correlations between the SA/V ratio of the printed objects and the mean dissolution time (MDT) of the printed dosage forms. Secondly, complete release profiles should be predicted based on correlations between the SA/V ratio and suitable mathematical equations. Validation experiments should be performed in all cases to appraise the quality of the prediction approaches.

Materials
Formulation 1, referred to as the PVA-PDM formulation, consisted of 5% (w/w) pramipexole dihydrochloride monohydrate (PDM, Chr. Olesen, Denmark) as an API of the biopharmaceutics classification system (BCS) class I, declared as good water solubility (c s ≥ 200 mg/mL) [49,50]. Mannitol (Parteck M ® , Merck, Germany) was used as a plasticizer at 10% (w/w) content. Polyvinyl alcohol (84%, PVA, Parteck MXP ® , Merck, Germany) was selected as a polymer. Formulation 2, referred to as the EVA-LD formulation, consisted of 10% (w/w) levodopa (Zhejiang Wild Wind Pharmaceutical, Dongyang, China), and was also a BCS class I API (c s ≥ 12 mg/mL) [51]. As water soluble component, a vinylpyrrolidone-vinyl acetate copolymer (VP-VA) was used (39.5%, Kollidon VA 64 ® , BASF, Ludwigshafen, Germany) and 10% mannitol was added as a plasticizer. The matrix consisted of ethylene vinyl acetate (EVA) with a content of 18% vinyl acetate (39.5%, Escorene ® FL 01418, TER Chemicals, Hamburg, Germany). To improve flowability, 1% fumed silica (Aerosil ® 200 VV Pharma, Evonik, Germany) was added to both formulations. Formulation 3, referred to as the PVA-PZQ formulation, consisted of 5% (w/w) praziquantel (PZQ, donated from Bayer AG, Leverkusen, Germany) as an API of BCS class II (c s = 0.4 mg/mL) [52][53][54], declared as poorly water-soluble. As a polymer basis, PVA was chosen with 95% content. All filament formulations were systematically developed to minimize diameter fluctuations of the filaments and to ensure highest printability with the available equipment. The criteria for the drug selection were heat stability and different classifications in the BCS. The melting point of PDM is also the decomposition point at 296-305 • C [20,[55][56][57]. LD melts and decomposes at 260-330 • C [58]. PZQ has its melting point already at 140-143 • C but decomposes only at temperatures >400 • C [59,60]. The investigated dose ranges do not correspond to therapeutic dosages and result from the drug loading of the filament and the volume of the objects. The polymer matrix of the first formulation should be water soluble and generate prolonged drug release. PVA fulfils both criteria as the polymer forms a water-soluble hydrocolloid matrix [61]. To test the transferability of the predictive model to other formulations, an inert, non-swelling matrix, EVA, was chosen [62]. To improve the printability and hydrophilicity of the filament, VP-VA was added. PVA was again chosen for transferring the model to the BCS class II drug, for the extended drug release as well as better printability. Thus, variations resulting from the process could be eliminated and the differences clearly attributed to the API.

Hot Melt Extrusion
The drug-containing filaments of formulation 1, 2 and 3 were produced via hot-melt extrusion with a co-rotating twin-screw extruder (Pharmalab HME 16, Thermo Fisher Scientific, Waltham, USA) using an in-house manufactured die with a diameter of 1.85 mm. The feed rate was set to 5 g/min and the screw speed to 30 rpm. The temperatures of the heating zones as well as the screw configurations are shown in Table 1. The filaments were hauled off with a winder (Model 846700, Brabender, Duisburg, Germany) at a speed of 1.8 m/min to the target diameter of 1.75 mm. The diameter was controlled using a laser-based diameter measurement module (Laser 2025 T, Sikora, Bremen, Germany).

3D Printing of Tablets
The drug-loaded filaments were printed on a FDM 3D printer (Prusa i3 Mk3, Prusa Research, Prague, Czech Republic) to oral dosage forms in various geometries. The geometries were designed with Fusion 360 ® (Autodesk, San Rafael, CA, USA) and sliced in Simplify3D ® (Simplify3D, Cincinnati, OH, USA) to obtain the desired G-code. The print temperature for the PVA-PDM filaments was set to 185 • C, the bed temperature to 60 • C and the printing speed was 20 mm/s. The printing temperature had to be increased for the EVA-LD formulation, as the filament was very flexible and could not be transported reliably through the nozzle by the conveying wheels in the print head at lower temperatures, as the filament would otherwise wrap around the wheels. To ensure a constant filament-flow through the nozzle, the temperature of the nozzle was set to 220 • C and the printing speed was reduced to 10 mm/s. For the PVA-PZQ formulation, a nozzle temperature of 185 • C could again be used, but the bed temperature had to be increased to 90 • C because the printed objects adhered poorly to the print bed and detached from the bed more often. Therefore, a printing speed of 10 mm/s was also selected here. To obtain high print accuracy, the layer height was set to 0.1 mm and the extrusion width to 0.4 mm using a nozzle with a diameter of 0.4 mm. The infill percentage of the concentric infill was set to 100%.  [63,64], release studies (n ≥ 3) were performed with the basket method (Method 1) in a dissolution tester (DT 726, Erweka, Langen, Germany). The baskets were 3D printed from water insoluble polylactide acid (PLA). They had to be adapted for printed tablets, since the mesh size of the regular Ph. Eur. baskets is small (0.36-0.44 mm) and the baskets were clogged by the swollen polymer of the PVA formulation. This affected the hydrodynamic medium flow around the printed oral dosage forms. The self-printed baskets have the same outer dimensions as the Ph. Eur. baskets, except that the mesh width was changed to 3 mm. The dissolution test for PDM containing tablets was performed in 500 mL of degassed 0.1 N hydrochloric acid at pH 1.2 under sink conditions (c s ≥ 200 mg/mL [50]; maximum concentration 0.08 mg/mL) and stirred at 50 rpm at a temperature of 37 ± 0.5 • C. The released API was measured using an UV-Vis spectral photometer (UV-1800 Shimadzu, Japan) at a wavelength of 263 nm. Dissolution testing of the levodopa containing geometries was performed in 1000 mL of degassed 0.1 N hydrochloric acid at pH 1.2 under sink conditions (c s ≥ 12 mg/mL [51]; maximum concentration 0.05 mg/mL) and stirred at 50 rpm at a temperature of 37 ± 0.5 • C. The API release was recorded with the same UV-Vis spectral photometer at a wavelength of 280 nm. The dissolution tests with PZQ were performed in 750 mL of degassed 0.1 N hydrochloric acid at pH 1.2 for the first 120 min and then transferred to 1000 mL phosphate buffer pH 6.8 with 250 mL of degassed 0.2 N tri-sodium phosphate dodecahydrate solution. The temperature was adjusted to 37 ± 0.5 • C. The drug release was measured at a wavelength of 210 nm. The tests were performed under sink conditions (c s = 0.4 mg/mL [54]; maximum concentration: 0.03 mg/mL). Samples were taken every 5 min for the first 30 min, then every 10 min for the next 90 min, followed by sampling in 20 min intervals. After a release time of 240 min, samples were taken every 30 min. The mean dissolution time (MDT) is used as a characteristic value that describes the drug release and was calculated according to Equation (1). It is expressed in units of time, usually in minutes, and indicates how long an API molecule is retained in a dosage form on average during dissolution [27,65].
ABC stands for the area between the curves and is calculated via the trapezoidal equation with c as the concentration of the API released over time t and c ∞ as the initial drug load of the dosage form. Values up to 100% of the release curve were used for the calculation since the ABC does not change afterwards.
Various mathematical models were tested to fit the complete release profiles for prediction. First, the Korsmeyer-Peppas equation, or power law, was applied (Equation (2)) [27,66].
The released amount of the API over time t is described by the term M t . The total amount of the API in the dosage form at the beginning of the release process is represented by M ∞ . The constant k describes the structural and geometrical characteristics of the dosage form (cylinder, sphere, film), also known as the reaction rate constant or release velocity constant. The diffusional constant n describes the underlying drug release mechanism. Only the first 60% of the release curve should be used for the calculation. The model can be used to analyse the underlying release properties of the system, when the mechanism is unknown, or more than one release mechanism is involved. It should be taken into account that the equation requires some properties of the matrix to be valid: the drug in the matrix must be distributed homogeneously, show unidirectional diffusion, and during dissolution, sink conditions should be maintained. Depending on the value of n, the release behaviour can be categorized into Fickian Diffusion, anomalous transport, or case-II transport (see Table 2) [29,66]. The second equation used for the description of the release profile was the Peppas Sahlin equation (Equation (3)) [29].
Using this equation, it is possible to describe the anomalous drug release process: the first term describes the Fickian diffusion and the second the Case-II relaxational contribution. The exponent n is, as in the Korsmeyer-Peppas equation, the diffusion exponent for any geometrical shape. The constants k 1 , k 2 describe the kinetics. The equation refers only to the anomalous release of the API [67].
The drug release from an inert matrix is often described using the Higuchi equation (Equation (4)). The Higuchi equation describes the release of a suspended drug from an insoluble matrix via diffusion and is also called square-root-of-time kinetic. This equation was initially intended for ointments but can also be used for solid dosage forms [68,69].
This calculation is used to calculate Q, the amount of API released at time t per unit area. It is composed of D, the diffusion coefficient of the API in the matrix, c 0 , the initial dose in the dosage form and c s , the saturation concentration of the API. The formula considers that the diffusion distance covered by the API does not remain constant but increases steadily during dissolution. The API concentration in the matrix decreases towards the medium since this is where the API is released first. Due to the growth of the diffusion distance, the concentration gradient decreases and, accordingly, the release rate. The Higuchi equation can be applied to suspension ointments and planar matrices (films) in which the API is dispersed as the equation describes a one-dimensional diffusion under sink conditions. The matrix should be inert and not undergo any change in structure, not even swelling, so that the diffusion coefficient remains constant. In addition, the drug concentration in the dosage form must be higher than the saturation solubility c s . As soon as these conditions are no longer given, the equation is invalid, since the API is then no longer dispersed but dissolved. The application of this equation is limited to 60% of the released API [70].
In this equation, M t represents the released API at time t and M 0 indicates the initial amount of the API in the tablet. The formula directly includes the SA/V ratio and the diffusion coefficient D of the API in the matrix.
For eroding tablets, the release can be described by the Hixson Crowell equation (Equation (6)). The formula reflects the size decrease of the dosage form. It assumes that the surface of the form decreases proportionally over time and the geometrical form remains constant. 3 W 0 reflects the initial amount of API and W i is the remaining API content over time t. The constant k is the constant of incorporation, which describes the surface area and volume relation of the dosage form. For the Hixson Crowell model, the release described is limited by dissolution velocity and not by diffusion of the API through the polymer matrix. Throughout dissolution, it is assumed that the surface decreases as layers detach over time [27,73,74].
Furthermore, a release of an eroding matrix can be described by the Hopfenberg equation (Equation (7)). The Hopfenberg equation is used to describe the drug release from different geometries (plates, spheres, and cylinders). The assumption is that drug release occurs via erosion and is not affected by diffusion.
The term M t M ∞ describes the release of the initial concentration c o of the system over time t with M t as API amount at the observed point of time and M ∞ as the total amount of API solute of the dissolved tablet at infinite time. The erosion constant is described by k 0 . The variable a 0 presents the radius of the selected geometry (sphere, cylinder, plate). The factor n represents the shape of the dosage form: n = 1 for a plate, n = 2 for a cylinder and n = 3 for a sphere. The limiting factor of this model is erosion of the matrix, which is the rate-limiting step of drug release and assumes that time-dependent diffusion resistances inside or outside the eroding matrix have no effect on it [27,75,76].
The Weibull equation (Equation (8)) can be used to describe release curve profiles, regardless of the underlying physical mechanisms (diffusion or erosion).
The value of c t is the concentration of the API at timepoint t. T i represents the lagtime before the release of API starts. If there is no delayed release, T i becomes zero. The parameter b describes the shape of the curve (exponential b = 1; sigmoid b > 1 or parabolic b < 1). The scale parameter a characterizes the time scale of the dissolution. As this equation is an empirical model without any kinetic basis, it does not provide any information about the underlying release properties, and only describes the curve profile. [27,[77][78][79].

Prediction
In order to be able to evaluate how accurate the prediction of the release curves is, the root mean square error of prediction (RMSEP) was calculated (Equation (9)) [42,80].
The RMSEP value describes the differences between the predicted and the measured values. The factor y i represents the experimental values andŷ represents the prediction values for every point of time of the release curve. The number of samples is indicated by n.

Comparison of the Dissolution Profiles
To compare the release curves, the similarity factor was used. This is approved by the FDA. The calculation was performed with Equation (10) [27,81,82].
In this equation, R t and T t stand for the mean released amounts of the API in % at time point t of the reference and the test product and n for the number of time points. A minimum of 12 measurement points was used as the mean to determine the f 2 value. The f 2 value is sensitive to the number of measurement points. For this reason, only one measurement point was included in the calculation after 85% released API. A f 2 value around 100 is desired, which indicates that the curves are identical. A value of 50 or more is accepted, which indicates that the values differ by max. 10%. Values below 50 indicate that the curves can no longer be considered similar.

Characterization of the Printed Tablets
To evaluate the precision of the printing and to see whether the designed SA/V ratio was also adhered to in the printed geometries, the printed tablets were measured with a caliper (CD-15CPXR, Mitutoyo, Kawasaki, Japan). The weight of the tablets was determined on a balance (Type 1702, Sartorius, Göttingen, Germany), and the tensile strength of the tablets was determined (n = 10) with the crushing force tester (TBH210, Erweka, Langen, Germany).

Characterization of the Printed Tablets
To ensure the printed tablets corresponded to the designed SA/V ratio, they were characterized (Supplementary Material, Tables S1-S7). The printing precision is considered very good for all three formulations and the deviations within a batch are very low, which indicates a good reproducibility. The largest weight variations were found for the PZQ-PVA formulation. The larger and heavier the printed tablet, the higher the variations in weight. In contrast, the API content is very similar, so the variations are small. These data suggest a good extrusion process with a homogeneous API distribution, as well as a robust printing process with small variations in precision. This leads to small variations in mass, which in turn leads to small variations in content. The printed PVA tablets have a high internal hardness, which means that they cannot be broken in the crushing force tester. The EVA tablets are very flexible, so that they also do not break in the crushing force tester but deform.

Drug Release from Dosage Forms with Defined SA/V Ratios
For initial release studies, various geometries of constant SA/V ratios with different doses were printed with the PVA-PDM formulation (Figure 1). The SA/V ratios were 1, 1.5, and 2 mm −1 , based on changes of the dimensions (Table 3). By keeping the SA/V ratio constant, the dose was varied according to the volume of the objects.

Characterization of the Printed Tablets
To ensure the printed tablets corresponded to the designed SA/V ratio, they were characterized (Supplementary Material, Tables S1-S7). The printing precision is considered very good for all three formulations and the deviations within a batch are very low, which indicates a good reproducibility. The largest weight variations were found for the PZQ-PVA formulation. The larger and heavier the printed tablet, the higher the variations in weight. In contrast, the API content is very similar, so the variations are small. These data suggest a good extrusion process with a homogeneous API distribution, as well as a robust printing process with small variations in precision. This leads to small variations in mass, which in turn leads to small variations in content. The printed PVA tablets have a high internal hardness, which means that they cannot be broken in the crushing force tester. The EVA tablets are very flexible, so that they also do not break in the crushing force tester but deform.

Drug Release from Dosage Forms with Defined SA/V Ratios
For initial release studies, various geometries of constant SA/V ratios with different doses were printed with the PVA-PDM formulation (Figure 1). The SA/V ratios were 1, 1.5, and 2 mm −1 , based on changes of the dimensions (Table 3). By keeping the SA/V ratio constant, the dose was varied according to the volume of the objects. The dissolution data ( Figure 2) revealed that dose and size did not influence the relative drug release but only the SA/V ratio, which was expected. The PVA matrix forms a hydrocolloid structure that swells first and slowly dissolves over time. The dissolution data ( Figure 2) revealed that dose and size did not influence the relative drug release but only the SA/V ratio, which was expected. The PVA matrix forms a hydrocolloid structure that swells first and slowly dissolves over time.  The smaller the SA/V ratio was, the longer it took until 100% of the API was released as the diffusion pathways are longer and vice versa. The geometries with a SA/V ratio of 1 mm −1 released 80% of the API in 100 min, the geometries with a SA/V ratio of 1.5 mm −1 released 80% in 60 min and with a SA/V ratio of 2 mm −1 , 80% API was released in 45 min. The dosages varied by a factor of 7.5 for the SA/V ratio of 2 mm −1 , by a factor of 4.5 for the SA/V ratio of 1.5 mm −1 and by a factor of 3.4 for the ratio of 1 mm − 1 .
The respective MDTs were calculated using Equation (1) and are listed in Table 3. Geometries with the same SA/V ratio have very similar MDTs that match the observed trend of the dissolution data. The MDT of geometries printed with a SA/V ratio of 1 mm −1 was 61 ± 3 min, with a SA/V ratio of 1.5 mm −1 35 ± 2 min and with a SA/V ratio of 2 mm −1 25 ± 1 min across all designs. The similarity of the curves was evaluated with the similarity factor (Equation (10)). The cylinder was used as a reference for each SA/V group. The f 2 -values are above 50 in each group and show that the curves are similar to each other.

Correlation between MDT and SA/V Ratio
Based on this data, a correlation between the SA/V ratio and the MDT was hypothesized. For a better understanding of the correlation between these two features, geometries with additional SA/V ratios were printed and the API release tested with the PVA-PDM formulation. For this, a SA/V ratio of 0.8 mm −1 was chosen as minimum and 6.0 mm −1 as maximum. The obtained data is shown in Figure 3a.

Correlation between MDT and SA/V Ratio
Based on this data, a correlation between the SA/V ratio and the MDT was hypothesized. For a better understanding of the correlation between these two features, geometries with additional SA/V ratios were printed and the API release tested with the PVA-PDM formulation. For this, a SA/V ratio of 0.8 mm −1 was chosen as minimum and 6.0 mm −1 as maximum. The obtained data is shown in Figure 3a. The linear equation of the regression was used to predict the MDT of 3D printed tablets with different SA/V ratios. The predictions were made exemplarily for four different SA/V ratios not used for the regression model and compared with the experimentally determined MDT of the corresponding printed geometries (Table 4). The predictions are very similar to the measured MDTs. The deviations are less than 5 min, which is acceptable with respect to small variations in the measured release curves. With this assessment, it is possible to categorize the SA/V ratios in terms of their release profiles.
Similar experiments were conducted for the EVA-LD formulation to assess the validity of the approach. Geometries with variable dosages and SA/V ratios of 0.9-6.0 mm −1 were printed, and the MDT was calculated from the obtained dissolution data (Figure 4). ), which could be linearized by log-transformation of both axes (R 2 : 0.9977). The linear equation of the regression was used to predict the MDT of 3D printed tablets with different SA/V ratios. The predictions were made exemplarily for four different SA/V ratios not used for the regression model and compared with the experimentally determined MDT of the corresponding printed geometries (Table 4). The predictions are very similar to the measured MDTs. The deviations are less than 5 min, which is acceptable with respect to small variations in the measured release curves. With this assessment, it is possible to categorize the SA/V ratios in terms of their release profiles.
Similar experiments were conducted for the EVA-LD formulation to assess the validity of the approach. Geometries with variable dosages and SA/V ratios of 0.9-6.0 mm −1 were printed, and the MDT was calculated from the obtained dissolution data (Figure 4).
A similar relationship was obtained as described earlier for the PVA-PDM formu- shown in Table 5. Again, predicted MDTs describe the obtained data well. The variations between the predicted values to the experimentally determined values are similar to the variations from the PVA-PDM formulation.
The approach for predicting MDT was also tested for BCS class II using the PVA-PZQ formulation ( Figure 5).    Additional geometries with defined SA/V ratios were printed and the MDTs resulting from the release curves were calculated. The measured and predicted MDTs are shown in Table 5. Again, predicted MDTs describe the obtained data well. The variations between the predicted values to the experimentally determined values are similar to the variations from the PVA-PDM formulation.
The approach for predicting MDT was also tested for BCS class II using the PVA-PZQ formulation ( Figure 5). Additional geometries with defined SA/V ratios were printed (not used for the regression model) and the MDTs resulting from the release curves were calculated. The measured and predicted MDTs are shown in Table 6. Additional geometries with defined SA/V ratios were printed (not used for the regression model) and the MDTs resulting from the release curves were calculated. The measured and predicted MDTs are shown in Table 6. Again, predicted MDTs describe the obtained data precisely. The variations between the predicted values and the experimentally determined values are similar to the variations from the PVA-PDM or EVA-LD formulation.
In the inert EVA matrix, which does not swell nor erode, the API can only be dissolved and transported by diffusion through the inert matrix into the dissolution medium. The API release is more prolonged compared to the PVA matrix. Considering geometries with a SA/V ratio of 1.5 mm −1 , the drug was released from the PVA-PDM matrix with a MDT of 35.58 ± 1.87 min, from the PVA-PZQ matrix with a MDT of 45.49 ± 3.27 min and from the EVA-LD matrix with a MDT of 111.07 ± 18.37 min. Despite the different behaviour of the matrix systems, the solubility properties of the APIs and the resulting differences in MDT, a correlation of the SA/V ratio and the MDT can be established for all three tested formulations. With the help of this correlation, the release time of an existing geometry with a given SA/V ratio can be predicted without having to waste material and perform dissolution tests. Similarly, the required SA/V ratio can be predicted for a desired MDT and the corresponding geometry can be designed and printed on this basis.
If it can be demonstrated that a linear relationship exists between ln(SA/V) and ln(MDT), a simple regression model can be generated using a two-point calibration with dissolution data only from geometries with a low and a high SA/V ratio. This would save further time and production costs and valid predictions could still be made.

Modeling and Prediction of Release Profiles
To predict the complete release curve and not only the MDT, different mathematical models were used to fit the data (Equations (2)-(8)) and the R 2 values of the fits were compared to evaluate the fit quality (Table 7). According to the Korsmeyer-Peppas model (KMP), the drug release from the PVA formulations follows an anomalous transport (diffusion constant n = 0.7 (PDM), n = 0.8 (PZQ), Table 2). Since anomalous transport can be described with different models (Equations (3) and (5)-(8)), only these equations were used to fit the data of the PVA formulations. For the EVA-LD formulation, n = 0.55 was obtained. This indicates a release according to square-root-of-t kinetic. This profile is often described with the Higuchi equation (Equation (4)). Additionally, the generally applicable Weibull equation (Equation (8)) as well as the Peppas Sahlin equation (Equation (3)) were used to fit the release profiles of the EVA formulation. For both BCS class I API formulations, the Peppas Sahlin equation provides the best fit. In order to also include enough data points for fast release dosage forms, the curves were fitted up to an API release of 98%. For the formulation with BCS class II API, the Weibull function provides the best fit ( Figure 6).
In earlier works, the slope of the Higuchi plot, the rate constant of dissolution and the constant k of the Korsmeyer-Peppas equation were already related to the particle size and the SA/V ratio of tablets [22,83,84]. In this study, all dissolution curves of BCS class I APIs were fitted with the Peppas Sahlin equation, as it provided the best fits and the constants k 1 and k 2 were determined (Table 8). For the PVA-PDM formulation, an average diffusion exponent of n = 0.79 was obtained, and for the EVA-LD formulation, an average diffusion exponent of n = 0.66 was obtained. The constants k 1 and k 2 were plotted in a graph as a function of the SA/V ratio. A linear correlation between k and the SA/V ratio was obtained by a log-transformation of both axes (Figure 7).
For both BCS class I API formulations, the Peppas Sahlin equation provides the best fit. In order to also include enough data points for fast release dosage forms, the curves were fitted up to an API release of 98%. For the formulation with BCS class II API, the Weibull function provides the best fit ( Figure 6). In earlier works, the slope of the Higuchi plot, the rate constant of dissolution and the constant k of the Korsmeyer-Peppas equation were already related to the particle size and the SA/V ratio of tablets [22,83,84]. In this study, all dissolution curves of BCS class I APIs were fitted with the Peppas Sahlin equation, as it provided the best fits and the constants k1 and k2 were determined (Table 8). For the PVA-PDM formulation, an average diffusion exponent of n = 0.79 was obtained, and for the EVA-LD formulation, an average diffusion exponent of n = 0.66 was obtained. The constants k1 and k2 were plotted in a graph as a function of the SA/V ratio. A linear correlation between k and the SA/V ratio was obtained by a log-transformation of both axes (Figure 7).    Using this relation and the resulting linear regression, k 1 and k 2 can be determined for arbitrarily selected SA/V ratios. If these calculated values are inserted into the Peppas Sahlin equation with the diffusion exponent n = 0.79, a concentration curve of the released API can be calculated for the specific time points. These calculations were performed exemplarily for the four different SA/V ratios not used for the model creation. The graphs comparing the experimental data points with the prediction curve are shown in Figure 8. Using this relation and the resulting linear regression, k1 and k2 can be determined for arbitrarily selected SA/V ratios. If these calculated values are inserted into the Peppas Sahlin equation with the diffusion exponent n = 0.79, a concentration curve of the released API can be calculated for the specific time points. These calculations were performed exemplarily for the four different SA/V ratios not used for the model creation. The graphs comparing the experimental data points with the prediction curve are shown in Figure 8.  The RMSEPs were calculated (Equation (9)) for the predictions and experimental results. The best predictive power was observed for a SA/V ratio of 1.6 mm −1 with a RMSEP of 0.56% and 0.9 mm −1 with a RMSEP of 0.74%. For the SA/V ratio of 2.3 mm −1 , the prediction above 70% API release underestimated the experimental data points, which results in a RMSEP of 2.85%. If only the release data points up to 70% released API are included, the RMSEP was 1.70%. Predicting the SA/V ratio of 4.7 mm −1 resulted in a RMSEP of 3.41%. If the prediction is only made up to 70% API release, the RMSEP is 1.18%, but only three time points are compared due to the quick drug release. Apparently, the predictions represent the experimental values better for smaller SA/V ratios because multiple time points can be considered and therefore the curves are not as sensitive to changes in k values. Nevertheless, the prediction approximated the experimental release curves well, especially in the early phase of drug release.
The prediction approach was tested for the EVA-LD formulation, respectively. Again, the constants k 1 and k 2 were plotted in a graph as a function of the SA/V ratio. A linear correlation between k and the SA/V ratio was obtained by a log-transformation of both values (Figure 9). Using this relation and the resulting linear regression, k 1 and k 2 can be determined for arbitrarily selected SA/V ratios. These calculations were exemplarily for the three different SA/V ratios, the MDT was already predicted in Section 3.2. The resulting graphs comparing the experimental data points with the prediction curve are shown in Figure 10.
Again, the constants k1 and k2 were plotted in a graph as a function of the SA/V ratio. A linear correlation between k and the SA/V ratio was obtained by a log-transformation of both values (Figure 9). Using this relation and the resulting linear regression, k1 and k2 can be determined for arbitrarily selected SA/V ratios. These calculations were exemplarily for the three different SA/V ratios, the MDT was already predicted in Section 3.2. The resulting graphs comparing the experimental data points with the prediction curve are shown in Figure 10.   The prediction approach was tested for the EVA-LD formulation, respectively. Again, the constants k1 and k2 were plotted in a graph as a function of the SA/V ratio. A linear correlation between k and the SA/V ratio was obtained by a log-transformation of both values (Figure 9). Using this relation and the resulting linear regression, k1 and k2 can be determined for arbitrarily selected SA/V ratios. These calculations were exemplarily for the three different SA/V ratios, the MDT was already predicted in Section 3.2. The resulting graphs comparing the experimental data points with the prediction curve are shown in Figure 10.   For the SA/V ratio of 1.89 mm −1 , the RMSEP value was 1.60% and for the SA/V ratio of 1.73 mm −1 , the RMSEP value was 1.95%. The RMSEP for the SA/V ratio 4.67 mm −1 results in a value of 3.42%. Again, the experimental data of the smaller SA//V ratios can be predicted better by the model, nevertheless, the prediction of higher SA/V ratios is still close to experimental data. The data implies that the developed prediction model also works for this formulation with an inert matrix.
For the PVA-PZQ formulation, the Weibull function provides the best fit (Table 7). Therefore, the approach for predicting the dissolution curves was tested for this equation. For the given SA/V ratio, the variables a and b were calculated (Table 9) and plotted with log-transformed axes ( Figure 11). The resulting linear equations were used to calculate the required variables for the SA/V ratios, whose curves were to be predicted. The comparison of the predicted and experimentally determined release profiles is shown in Figure 12.  The sigmoidal curve is described with the Weibull function. Initially, the amount of API released is somewhat underestimated, but as dissolution progresses, the mean value is predicted well. For the curves of lower SA/V ratios, the course of the release curve is overestimated from 55-80% released API and a stronger S-curve course is predicted. With increasing SA/V ratio, however, the Weibull curve predicts the release profile better. The RMSEP for the SA/V ratio of 1.3 mm −1 is 2.5%, for 1.83 mm −1 3.6%, for 2.3 mm −1 2.1% and for 4.67 mm −1 1.0%. The predicted value of the curve deviates from the experiment by 2.3% on average, which can be considered acceptable.
With this approach, it is possible to extend the prediction of the MDT as described in  The sigmoidal curve is described with the Weibull function. Initially, the amount of API released is somewhat underestimated, but as dissolution progresses, the mean value is predicted well. For the curves of lower SA/V ratios, the course of the release curve is overestimated from 55-80% released API and a stronger S-curve course is predicted. With increasing SA/V ratio, however, the Weibull curve predicts the release profile better. The RMSEP for the SA/V ratio of 1.3 mm −1 is 2.5%, for 1.83 mm −1 3.6%, for 2.3 mm −1 2.1% and for 4.67 mm −1 1.0%. The predicted value of the curve deviates from the experiment by 2.3% on average, which can be considered acceptable.
With this approach, it is possible to extend the prediction of the MDT as described in Section 3.2 and predict the release profile of different APIs of various BCS classes for a given SA/V ratio. As shown with the different formulations and APIs, this method can potentially be applied for various matrices and drug substances. This approach enables more accurate predictions to be made for the release behaviour of 3D printed dosage forms.

Discussion
A correlation between MDT and SA/V ratio was found in these experiments and used for a prediction of drug release profiles with three different formulations. With the help of this approach, the MDT of printed dosage forms with known SA/V ratios can be easily predicted and it is possible to categorize the SA/V ratios in terms of their release. MDT predictions are quite accurate for all three formulations, varying from experiment mostly by ≤5 min, which is acceptable. The largest variations are found for the EVA formulation, which may be due to the more difficult printing conditions. Due to the high flexibility of the filaments, the constant transport of the filament through the nozzle is not always given, resulting in less accurate printing of the tablets. In addition, the release of the API is longer than in the PVA formulations, resulting in larger standard deviations. Since inert matrices often lead to prolonged release kinetics and, therefore, require time consuming experiments, the identified relationship between MDT and SA/V ratio represents a helpful tool for dosage form design. Accurate predictions of release profiles of BCS class I APIs could be made using the Peppas Sahlin equation. This equation was used to describe the dissolution curves of two different matrix systems with different release characteristics. The constants of the equation were found to change in a predictable manner with the SA/V ratio. For both formulations, precise predictions of release profile could be made, with an average variation of 0.6-3% of API released per time point. Predictions for smaller SA/V were mostly more accurate than those for larger SA/V, which may be due to the number of time points included, making the curves more sensitive to changes in the constants. In the PVA formulation with the BCS class II compound, the Weibull function was found to be most suitable to describe the drug release (R 2 : 0.999). Again, the linear relationship between the logarithmised constants and the SA/V ratio was used to calculate dissolution predictions. With an average deviation of 2% from the predicted to the experimentally determined release profile, it can be stated that dissolution curves of dosage forms with a poorly water-soluble API can also be modeled.
Future studies should include testing more formulations and drug substances to investigate the validity of this approach for other drugs and polymers, and expanding the database, as well as investigating the use of ANNs for these kinds of dosage forms. In addition, it should be tested whether the prediction model is also transferable to tablets or dosage forms manufactured by other 3D printing processes, e.g., by powder bed 3D printing or semi-solid 3D printing with different materials and porosity.

Conclusions
We demonstrated that the relationships between constants of established equations and the SA/V ratio can be exploited to predict complete dissolution curves for specific SA/V ratios. The validity of this approach was shown for water-soluble as well as poorly water-soluble drug substances and erodible and inert polymers. It might be possible to establish a dataset for polymer matrices with various APIs of different BCS classes and directly calculate the release curves of dosage forms with a known SA/V ratio. Such an approach would save additional costs and time-consuming dissolution studies of individual geometries and matrix systems, which is especially important for community pharmacies and hospitals. Future patient care with 3D printed oral dosage forms could be: 1. the optimal drug release for a patient is determined, 2. the appropriate formulation and SA/V ratio are selected, 3. a suitable geometry is 3D printed with the required dose of the API. Therefore, with an established matrix system, not every drug release of oral dosage forms would need to be tested in dissolution studies but could be mathematically predicted. As a result, the dissolution rate could be individualized in decentralized settings, an approach that has not been easily adaptable before.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/pharmaceutics13091453/s1, Table S1: Physical characterization of the printed tablets for SA/V ratio 1-2 mm −1 (n ≥ 3, x ± s), Table S2: Physical characterization of the printed tablets for the correlation generation with the PDM-PVA formulation (n ≥ 3, x ± s), Table S3: Physical characterization of the printed tablets for the correlation generation with the LD-EVA formulation (n ≥ 3, x ± s), Table S4: Physical characterization of the printed tablets for the correlation generation with the PZQ-PVA formulation (n ≥ 3, x ± s), Table S5: Physical characterization of the printed tablets for the prediction validation with the PDM-PVA formulation (n ≥ 3, x ± s), Table S6: Physical characterization of the printed tablets for the prediction validation with the LD-EVA formulation (n ≥ 3, x ± s), Table S7: Physical characterization of the printed tablets for the prediction validation with the PZQ-PVA formulation (n ≥ 3, x ± s).