Optimizing the Entrainment Geometry of a Dry Powder Inhaler

12 Pages·2016·1.75 MB·English

by  Thomas Kopsch

PharmRes
DOI10.1007/s11095-016-1992-3
RESEARCHPAPER
Optimizing the Entrainment Geometry of a Dry Powder Inhaler:
Methodology and Preliminary Results
ThomasKopsch1 &DarraghMurnane2&DigbySymons1
Received:25April2016/Accepted:30June2016
#TheAuthor(s)2016.ThisarticleispublishedwithopenaccessatSpringerLink.com
ABSTRACT KEY WORDS boundary-condition .cost-function .DPI .
Purpose Forpassivedrypowderinhalers(DPIs)entrainment entrainment .optimization
andemissionoftheaerosolizeddrugdosedependsstronglyon
devicegeometryandthepatient’sinhalationmanoeuvre.We
proposeacomputationalmethodforoptimizingtheentrain- ABBREVIATIONS
mentpartofa DPI. The approach assumes that thepulmo- BC BoundaryCondition
nary delivery location of aerosol can be determined by the CFD ComputationalFluidDynamics
timingofdoseemissionintothetidalairstream. COPD ChronicObstructivePulmonaryDisease
Methods An optimization algorithm was used to iteratively DPI DryPowderInhaler
perform computational fluid dynamic (CFD) simulations of EE Eulerian-Eulerian
thedrugemissionofaDPI.Thealgorithmseekstoimprove EL Eulerian–Lagrangian
performance by changing the device geometry. Objectives KTGF KineticTheoryforGranularFlow
were to achieve drug emission that was: A) independent of PIF PeakInhalationFlow
inhalationmanoeuvre;B)similartoatargetprofile.Thesim-
ulationsusedcompleteinhalationflow-rateprofilesgenerated
dependent on the device resistance. The CFD solver was INTRODUCTION
OpenFOAM with drug/air flow simulated by the Eulerian-
Eulerianmethod. EstablishmentofDesignObjectivesforInhaledDrug
Results Todemonstratethemethod,a2Dgeometrywasop- Delivery
timized for inhalation independence (comparing two breath
profiles) and an early-bolus delivery. Entrainment was both Dry powder inhaler (DPI) therapy has increased in impor-
shear-drivenandgas-assisted.Optimizationforadelayinthe tance asa therapeuticmodality forbothpulmonary diseases
bolusdeliverywasnotpossiblewiththechosengeometry. such as asthma, and systemic diseases such as diabetes. In
Conclusions ComputationaloptimizationofaDPIgeometry particular, the adoption of the Montreal Protocol banning
for mostsimilar drugdeliveryhas beenaccomplished for an the use of chlorofluorocarbon propellants accelerated devel-
exampleentrainmentgeometry. opmentofDPIsasanalternativetothedifficultreformulation
ofpressurizedmetereddoseinhalers(1,2).Inrecentyears,the
suitabilityofDPIsforaerosolizationanddeliveryofhighdrug
payloads has been exploited, including the introduction of
inhaled insulin and antibiotic products (3,4). There has also
* ThomasKopsch been an interest in developing bioequivalent (generic) DPI
[email protected]
productsasmoreestablishedproductsreachtheendofpatent
exclusivity.Thereisthereforeaneedtodesignandoptimize
1 DepartmentofEngineering,UniversityofCambridge,Trumpington DPIdevicestoachieveavarietyoftherapeuticoutcomesina
Street,CambridgeCB21PZ,UK diverse patient population ranging from those with healthy
2 DepartmentofPharmacy,PharmacologyandPostgraduateMedicine, respiratory function to those with compromised respiratory
UniversityofHertfordshire,CollegeLane,HatfieldAL109AB,UK function.
Kopsch,MurnaneandSymons
WhendesigningaDPI,itisimportanttounderstandhow into the airstream correlated (16) with the region to which
particular design elements influence the release (emission) of drugwasdeposited.
aerosolizeddrugandconsequentlythedeviceperformance.It Consequently, in this study it is assumed that the lung
isthereforenecessarytopredictdrugreleaseprofilesforpar- location to which drug will be delivered correlates well
ticularentrainmentgeometriesandtorelatetheseprofilesto withthepointintheinhaledtidalvolumeintowhichdrug
thepharmaceuticalobjectivesofdrugdelivery(twoexamples is introduced during inhalation. Although not
include:achievingbioequivalence,ormaximizingperipheral disregarding the importance of the de-agglomeration pro-
airway deposition for systemic drug delivery). When cesses, the argument is valid in cases where the DPI de-
predictingdrugreleaseprofilestheentrainmentmechanisms signer cannot easily influence the size of drug particles, if
mustbeunderstood.Ingeneral,authorsdistinguishbetween we assume a powder that has been optimized for its ease
gas-assisted,shear-forcedriven,capillaryandmechanicalflu- ofdispersion(17,18). Followingthisargument,itisimpor-
idization(5).Thefluidizationmechanismsfordifferentdevices tanttooptimizeDPIstoreleasedrugatthecorrectinstant
maybefoundinliterature(e.g.(6,7)). with respect to the fraction of inhaled volume (x) rather
When designing an inhaler, several objectives have to be than with respect to the inhalation time (t).
considered: First, drug should deposit in similar pulmonary
airways, independent of the patient’s pulmonary function
and ability touse the device (objectiveA).For example,van ComputationalFluidDynamicsforDPIDesign
den Berge (8) stresses the importance of targeting the small
airwaysforeffectivetreatmentofasthmaandchronicobstruc- Computationalfluiddynamic(CFD)analysisreferstotheuse
tive pulmonary disease (COPD). Second, the drug release ofcomputationaltechniquestosolvefluidflowproblems.The
profile should be similar to a (pre-defined) desired released application of CFD techniques in inhaled drug delivery has
profile, for example an early-bolus release (objective B) (9). increasedinrecentyearsforpredictionofpulmonaryaerosol
Additionalobjectivesaretoachievede-agglomerationofdrug deposition, but also to improve understanding of aerosoliza-
powder (10,11) and low-manufacturing costs but these are tion processes, e.g. to predict the performance of DPIs (19).
beyondthescopeofthispaper. CFD has also been used to assess device resistance (20,21),
Itisacceptedthattheparticlesizedistributionoftheemit- track the path of individual particles (21) and to investigate
ted aerosol is a significant factor in determining the site of particlede-agglomerationinDPIs(10,22,23).
deposition in the pulmonary airways (12), and the degree of InCFDanalysisforDPIstheparticulate(drug)phasecan
de-agglomerationduringdeviceactuationiscrucialindeter- eitherbeignored(asinglephaseapproachconsideringairflow
miningtheparticlesize.However,theemission-timeprofileof only)ormodelledusingamultiphaseapproach.Themajority
that aerosol (even one with a size distribution suitable for of reported studies using multiphase approaches for simula-
targeting delivery to specific airways) is an area of DPI and tion of DPIs have used Eulerian–Lagrangian (EL, particle-
inhaler design that remains under-explored. With the tracking)approaches(24).ELapproachesmakeitdifficultto
AdaptiveAerosolDeliverySystem(13),peripheralairwayde- studythefluidizationprocessofapowderbedduetothelarge
positionofaerosolproducedusingavibratingmeshnebulizer number of particle-particle interactions. It is for this reason
canbeenhancedbylimitingemissionandinhalationofaero- thatwehaverecentlyintroducedEulerian-Eulerian(EE)ap-
solto50–80%oftheinspiratorycycle.Anumberofauthors proaches(9,25)forinvestigationofDPIs.
arguethatdeliveringdrugtoaparticularlocationofthelung IntheEEapproach, thedrugparticlesare modelled asa
correlateswelltothepharmacokineticprofile(14).Ithasbeen secondcontinuousphaseandtheinteractionbetweenthegas
shownthatthetimingofreleasecorrelateswiththedeposition andthegranularphaseismodelled.Anewvariable,thevol-
locationfornebulizers.Scheuchetal.(15)foundthatthesiteof umefraction(α)isintroducedtokeeptrackofthelocalfrac-
drugdepositionmaybeinfluencedbythetimingofthebolus tionofdrugphase.EEapproachesareoftenusedinthestudy
delivery. Heyder et al. (16) investigated in which ‘volumetric of other fluidization processes (26,27). As in other CFD ap-
lungdepth’abolusofdrugdeposits.Drugwasintroducedinto proachesthedomainofinterestissplitintoalargenumberof
theinhaledairflowafteraparticularfractionofthetotalvol- cells. For each cell the Eulerian-Eulerian solver solves for a
ume was inhaled. Heyder et al. (16) showed experimentally numberofvariablesincludingα,theaveragevelocityofdrug
that particle deposition was higher when the aerosol bolus particles v and the average velocity of the gaseous phase
drug
penetrateddeepintothelung.Similarly,Brandetal.(12)de- v . As shown in Fig. 1 the value of α indicates the number
air
termined in which ‘lung volume element’ drug is likely to ofparticlesinthatcell.IntheEEapproach,allphases,includ-
depositbydeterminingtheprobabilityofaparticledepositing ing theparticulate phases, are modelled as continuous fluids
in the lung as a function of volumetric lung depth for an andconsequentlymass,momentumandenergyconservation
aerosol with a constant drug concentration. In the case of a lawsareappliedtoallphases.Inmanycases,theseequations
bolus delivery the instant at which the drug was introduced look similar to the standard Navier–Stokes equations. The
OptimizingdryPowderInhalers:Methods+Results
Fig.1 Idealizedrepresentationof
theEEapproach.Thevolume
fractionαindicatesthedensityof
particles.
kinetictheoryforgranularflow(KTGF)(28)isonemethodto drug emission performance of a simple entrainment
model the constitutive behaviour of a granular phase in the geometry.
EE approach. KTGF makes predictions about stresses and
phase interaction terms in the governing equations of the
fluid. METHODS
For modelling purposes, the applied boundary condi-
tion at the outlet of the inhaler is usually a constant flow CostFunctionDevelopment
rate or a transient pressure boundary condition (25). A
constant flow rate boundary condition is prescribed in Thegoalofnumericaloptimizationistofindanargumentx
regulatory inhaler testing (29,30). However, physiological- that minimizes a cost function C(x), where x consists of a
lyrelevantinhalationflow-rateprofilesaremorecomplex. number of design variables. In engineering applications the
As shown by de Koning et al. (31,32) the peak inhalation cost C measures some quality of a design. For example, C
flow (PIF), total inhalation time and other inhalation pro- could measure the manufacturing cost or energy consump-
file variables depend on the device resistance andthe type tion.Inthisstudyweonlyconsiderdrugdeliveryperformance
of disease. It is of interest to develop a DPI modelling ofaDPIdevice.Weproposetwocostfunctionsthatmeasure
approach that can incorporate the complexity of inter- entrainment and emission performance. The optimization
individual variability in inhalation performance and goal is to minimize these functions. Cost function A prefers
flow-rate profiles. For this reason the experimental data devicesthatcanemitdrugtosimilarpulmonaryairways(12)
from de Koning et al. (31) was employed to model (i.e.volumetriclungdepth),independentofthepatient’srespi-
resistance-dependent inhalation profiles and to establish ratoryfunction.CostfunctionBprefersdevicesthatachievea
the outlet boundary conditions for a CFD simulation. In desired drug release profile (e.g. early bolus). The cost func-
particular it is of interest to simulate the drug release tions are developments of earlier concepts presented by the
profiles during the simulated inhalations in order to guide authorsofthisstudy(9,25).
the design and optimization of different DPI geometries. Forthefirstdesignobjective,thegoalistoachievethe‘most
There has been relatively little prior work on the use of similar drug delivery’, i.e. the drug should penetrate to the
computational optimization methods for the design of DPIs. samepulmonaryregionfordifferentpatients.Totestthisob-
Theaimofthisstudywastoaddresstheknowledgegapinthe jective a large number of patients with different inhalation
applicationofcomputationalfluiddynamics(CFD)inthenu- manoeuvres could be considered. However, to illustrate the
merical optimization of DPIs. In order to achieve the aim, approachonlytwoinhalationprofiles,1and2,areconsidered
optimizationcostfunctionshavebeendevelopedforthequan- inthisstudy.
titativeassessmentoftheperformanceofdeviceentrainment In general, patients have different inhalation profiles; i.e.
geometriesforarangeofsimulatedinhalationprofiles.Asan whentheyinhale through the device, they produce different
example ofthe methodology, these cost functions have been volumetricflowratesQ (t)andQ (t)(seeexampleprofilesin
1 2
appliedinanoptimizationframeworkinordertoimprovethe Fig. 2a) leading to different drug emission rates (dM/dt)
Kopsch,MurnaneandSymons
Fig.2 Idealizedrepresentationsofinhalationflowrate(Q),inhaledvolume(VðtÞ),releaseddrug(M)andscaledvolume(x)asfunctionsoftimeoreachother.
Sampledataforclarificationpurposesonly.ThisFig.issimilartoFig.1in(1).(a)Flowrateandreleaseratevs.time.(b)Volumeofinhalationvs.time.(c)
Releaseddrugvs.volume.(d)Releaseddrugvs.scaledvolume.
(Fig. 2a). The integral (Eq. 1) gives the inhaled volume as a The devicewillideallyhavereleasedthesame amountof
functionoftime,t: drugforpatients1and2whenthesamefractionxofthetotal
Z (cid:2) (cid:3) inhaledairV hasbeeninhaled.Itfollowsthatidenticaldose
t i,Tot
V ðtÞ¼ Q t0 dt0 ð1Þ emissionwouldberepresentedbyEq.4,forallx.
i i
0 (cid:4) (cid:5) (cid:4) (cid:5)
M1 xV1;Tot ¼M2 xV2;Tot ð4Þ
wherei=patient1or2.Thetotalinhaledvolume(seeFig.2b)
isgivenbyEq.2: Since it is impossible to achieve Eq. 4 exactly, the design
Z ∞ (cid:2) (cid:3) objectiveis,therefore,tominimizethedifferencebetweenthe
Vi;Tot ¼ Qi t0 dt0 ð2Þ drugreleaseprofiles(Eq.5)forallx:
0 (cid:6) (cid:4) (cid:5) (cid:4) (cid:5)(cid:6)
After carrying out a CFD simulation of entrainment, the cðxÞ¼ (cid:6)M1 xV1;Tot −M2 xV2;Tot (cid:6) ð5Þ
emitted dose profile M can be calculated. M is the mass of
i i
drugthathasleftthedeviceatonepointandcanbeexpressed Figure2dshowsthesamereleaseprofilesasinFig.2c,but
aseitherafunctionoftimeM (t)orinhaledvolumeM (V)(see as a function of scaled volume x. It is hoped to minimize
i i
Fig.2c). M is the total drug amount in the primed dosing c(x) = 0 for all x. A possible cost-function to achieve this is
Tot
unit(i.e.capsule,blister,meteringcup).M (V)maybescaled thereforeEq.6:
i
alongtheV-axisbydefiningthescaledvolume(Eq.3): Z Z
1 1(cid:6) (cid:4) (cid:5) (cid:4) (cid:5)(cid:6)
V CA ¼ cðxÞdx ¼ (cid:6)M1 xV1;Tot −M2 xV2;Tot (cid:6)dx ð6Þ
x ¼ i where 0≤x≤1 ð3Þ 0 0
Vi;Tot
MinimizingC willensure‘mostsimilardrugdelivery’for
A
ThisgivesM(V)=M(xV ),(seeFig.2d). different patients. Note that since CFD solves discretized
i i i,Tot
OptimizingdryPowderInhalers:Methods+Results
versions of the governing equations in practice the integral have not released all the drug dose before PIF is reached
(Eq.6)hastobeapproximatedasasummation. (34). A fixed simulation time is therefore only acceptable if
Minimizing the cost function C guarantees a device air- mostofthedrugleavesthedevicewithinthesimulatedtime.
A
flow geometry that can achieve the highest similarity of the Apossiblesolutiontotheseproblemsistospecifytheflowrate
emission-volume profile. However, it cannot be guaranteed Q (t) as a function of time at the outlet and model a whole
thatthisemissionprofileisalsoadesiredone.Itdependson inhalationprofile(oratleastthefirstpartoftheprofileuntil
thedesiredtherapeuticapplicationwhether,forexample,ei- drugemissioniscomplete).
ther an early bolus delivery (to reach peripheral airways) or The appliedQ (t) BCshould beasrealisticaspossible.In
continuousemission(toachievedepositionthroughoutallair- practice measured flow rate profiles vary significantly as a
ways)isrequired.Weshould,therefore,alsoconsidercompar- functionofdeviceresistanceR,i.e.alowerflowrateisexpected
ingthereleaseprofilewithadesired,idealizedreleaseprofile. tobeproducedbyanygivenpatientforhighresistancedevices
As an illustration, if for a particular medical application an andviceversa.Forexample,deKoningetal.(31,32)plottedthe
earlybolusdeliveryofdrugisdesired,anidealisedearlybolus flow rate as a function of time for two different device resis-
profile(seeFig.2d)canbedefinedusingEq.7: tances and health of patient. Several important (averaged)
( ) parameters of the inhalation profile were reported, such as
x
M ðxÞ¼ Mtot x if x < x0 ð7Þ PIF, the time when the PIF occurs tPIF, and total time of
3 0
Mtot if x≥x0 inhalation ttot as a function of resistance and patient health.
Thus, it is possible to fit a polynomial to these points. The
This is a ‘ramp’ function, i.e. M (x) represents an early- followingprocedurewasapplied(seeFig.3a):
3
bolus profile, which increases linearly from 0 to x and is
0 pffiffiffiffiffi
constantfromx to1.Thisismorerealisticcomparedto(9), 1. DeterminetheresistanceRinQR¼ Δpbyrunninga
0
whereanidealstepfunctionwasused.Toevaluatehowclose- steady-state CFD simulation (using simpleFOAM from
ly the emitted dose profile from our device matches the OpenFOAM (35)) where a constant pressure difference
idealised profile the cost function given by Eq. 8 has been Δp is applied. Q is found from the simulation. Hence, R
pffiffiffiffiffi
used: ¼ Δp=Q iscalculated.
Z 1 Z 1(cid:6) (cid:4) (cid:5) (cid:6) 2. InterpolatethedatafromdeKoningetal.(31,32)tofind
CB ¼ cðxÞdx ¼ (cid:6)M1 xV1;Tot −M3ðxÞ(cid:6)dx ð8Þ relevantflowfieldparameters(PIF,tPIF,ttot)forthepartic-
0 0 ularvalueofRandpatients’diseasestate.
3. ModelaflowratefunctionQ(t)fromtheseparameters.
Minimizing cost function C gives a geometry with the
B
desiredreleasebehaviour.Toalsoachieve‘mostsimilardrug 4. Accountforafractionofbypassingair:ManyDPIdevices
delivery’ between different patients requires simultaneous allowthemajorityofinhaledairflowtobypassthedrug
entrainmentgeometry,e.g.(21).
minimizationofcostfunctionC .
A
Inprinciple,awiderangeoffunctionscouldbefittedtothe
BoundaryConditions
parametersinstep3tomodelaflowratefunctionQ(t).Inthis
case,however,Q(t)wasmodelledusingtwopiecewise-defined
WhenconductingaCFDanalysis,boundaryconditions(BCs)
3rdorderpolynomialsf(t)andg(t):
have to be specified. At the inlet a (constant) atmospheric
8 9
pressureisusuallyspecified.Theappliedboundarycondition
<fðtÞ if t < t =
attheoutletshouldresembletheflow-fieldgeneratedbyareal PIF
QðtÞ¼ gðtÞ if t≥t ð9Þ
patientinhalingthroughthedevice. : 0 if t≥PtIF ;
tot
In many experimental and computational studies a con-
stant flow rate Q or a constant pressure drop Δp is applied
(10,21,33).Aconstantflowrateboundaryconditionsimplifies Thefirstpolynomial,f(t),linkstheonsetoftheinhalation
experimentalstudiesandregulatorytesting(29,30).However, with the PIF. The second polynomial connects the PIF with
if much of the drug entrainment occurs before the PIF is theendoftheinhalation.Thecoefficientsofthepolynomialsf
reached,orifthePIFisnotsustainedforverylong,thenusing (t) and g (t) can be determined numerically by imposing a
aconstantflowrateBCmaynotmodeltherealentrainment numberofconditionsatthestartandendoftheinhalation:
behaviour accurately. Zimarev et al. (25) went some way to
addressingthisissuebyapplyingapressurerampBCofcon- 1. Att=0:Q(0)=f(0)=0
stantdpforafixedsimulationtimetorepresentthestartofan 2. Att=t :i.e.Q(t )=g(t )=0
dt tot tot tot
inhalation.Thisisagoodimprovementifentrainmentoccurs 3. Q (t) drops slowly towards the end of the inhalation, i.e.
early in the inhalation manoeuvre. However, some devices Q′(ttot)=0.
Kopsch,MurnaneandSymons
Fig. 3 Determiningthe boundarycondition and examples ofrelease profiles for different resistance R. (a) Determine transient boundary conditions.(b)
ModelledQ(t)functionfordifferentresistancesR -R.PiecewisedefinedpolynomialswerefittedtodatafromdeKoningetal.(31,32).
3 8
At the transition of the two polynomials, i.e. at point was generated depending on these design parameters.
(t ,PIF),itisrequiredthat BlockMesh(OpenFOAM(35))wasusedtomeshthegeome-
PIF
try. The Eulerian-Eulerian solver twoPhaseEulerFoam
4. Q(t)iscontinuous,i.e.f(t )=g(t )=PIF (OpenFOAM)wasusedtosimulatethedrugreleasethrough
PIF PIF
5. Q(t )isamaximum,i.e.f′(t )=g′(t )=0. the device for high and low flow-rate boundary conditions.
PIF PIF PIF
6. Itwasalsoassumedthatthesecondderivativeiscontinu- Transient boundary conditions were applied using the
ous,i.e.f″(t )=g″(t ). OpenFOAM library swak4Foam (37). A Python-script (36)
PIF PIF
determinedthedrug-releaseprofilesasafunctionofthescaled
Usingtheseconditions,thecoefficientsofthepolynomials inhaledvolumexforeachboundaryconditionandthencal-
weredetermined.Fig.3bshowsexamplefitstodatafromde culatedthecostfunctionsC andC .
A B
Koningetal.(31,32)fordifferentdeviceresistances. A total cost function C = w C + w C was calculated
tot A A B B
In this study the method described above, using de (where w and w are weightings). The cost function C
A B tot
Koning’sdata(31,32),wasusedtogeneratetheidealizedin- wasminimizedbysystematicallychangingthegeometryusing
halation profile Q (t) for an average patient for each device asteepestgradient-descentmethod,seee.g.(38).
1
geometryconsidered.ThesecondinhalationprofileQ (t)was
2
generatedinasimilarway,butwithalowerPIF.ThePIFof DefinitionofDesignSpaceforOptimizationandCFDStudies
Q ðtÞwaschosentobeonestandarddeviationbelowthatof
2
Q (t).ThetimetoPIFt andthedurationt werethesame AsimpleDPIdevicewasconsideredinthisstudy(Fig.4).The
1 PIF tot
forbothQ ðtÞandQ ðtÞ(incontrasttotheexampleprofiles geometrywaschosensothatitcouldrepresentarangeoftypes
1 2
showninFig.2). of design of existing DPI. The device consists ofan inlet, an
entrainment area and an outlet. A fixed volume of drug of
3 mm2 was initially placed into the entrainment area (note
CFDandOptimization that in 2D simulations, a volume may be represented by an
areaA).Thesizeoftheoutletisa,theheightofthecapsule/
CFDMethodology blisterandtheinletish.Thecapsule/blistermayormaynot
be filled with drug completely. The separation between the
Asanexampleimplementation,thecostfunctions(section2.1) drugsurfaceandtheupperboundaryofthecapsuleisc.dis
and thedynamically generatedboundary conditions (section theanglebetweentheinlettubeandthecapsule.Theamount
2.2)wereappliedinaCFDoptimizationroutine.First,asim- ofdrugAfilledintothedevicewasconstant.Theheight,h,of
ple 2D geometry was parametrized in terms of five design the capsule was therefore calculated from
variables:a;b;c;d;e.UsingaPython-script(36)thegeometry h¼cþA=ðeþaþbÞ.
OptimizingdryPowderInhalers:Methods+Results
OptimizationAlgorithm
In this study a steepest gradient descent optimization algo-
rithm was used to minimize the total cost function C .
tot
Successfulimplementationofgradientdescentreliesondeter-
miningthesearchdirectionpwhichisthedirectionofsteepest
descent p=−∇C and on deciding a suitable step length α
tot
(38). To enhance the performance of the optimization algo-
rithm, all five design variables were scaled to confine values
between0and1.Forexample,a*,thescaledversionofdesign
Fig.4 Parameterizedentrainmentgeometry.
variablea,wascalculatedusingEq.10.
a−a
a* ¼ low ð10Þ
Table I shows the lower and upper boundaries of the a −a
high low
design variables. The values of these boundaries were
chosentogivegeometriesofrealisticsizes.DaniherandZhu alowandahigharethelowerandupperboundaryofa,respec-
(5)identifieddifferenttypesofdrugentrainmentinDPIs,i.e. tively. Similarly, all other design variables were scaled. This
shear-driven and gas-assisted entrainment. This parameter- givesavectorofscaleddesignvariables:
(cid:4) (cid:5)
izeddeviceisdesignedtoallowdifferententrainmentbehav-
x¼ a*;b*;c*;d*;e* ð11Þ
iourfordifferentvaluesofthedesignvariables.Forexample,a
gas-assisted entrainment may be expected when e and c are
If, at step k of an optimization process, the current set of
Å
smallanddiscloseto90.Suchageometryapproximatesthe design variables is xk=(a*,b*,c*,d*,e*), a new choice of
designvariables,forthenextstep(k+1),maythenbeevalu-
bathtubshapedblisterfoundinsomeDPIswherethefoillidis
atedasfollows(38):
pierced to fromaninletand outlet(39). Onthe other hand,
shear-driven entrainment, e.g. (40), may be expected when e xkþ1 ¼ xkþαp ð12Þ
andcarelarge.Inaddition,theresistancemaybeinfluenced
byvaryinganyoftheparameters,sincethelatterchangethe In this study, the search direction p was determined by
airflowpathswithintheentrainmentregionoftheDPIdevice. approximating −∇C using a finite difference method.
tot
Usingabacktrackinglinesearchasuitablesteplengthαwas
determined(38).Sincethisisaconstrainedoptimizationprob-
Eulerian-EulerianModelling lem a penalty term was applied if the constraints were
exceeded.Inordertoavoidconvergencetoalocalminimum,
Whenthevolumefractionofthedrugphaseishigh,thedrug the same optimization algorithm was repeated for different
may significantly influence the flow field. For this reason an (random)initialvaluesofthedesignvariables.
Eulerian-Eulerian (EE) approach was used in this study to
simulate drug release profiles. When the volume fraction of
drugisveryhigh,theEEapproachismore computationally
efficient than an Eulerian–Lagrangian (EL) particle tracking RESULTS&DISCUSSION
approach.Themodelsandparametersusedinthisstudyare
ExperimentalValidationoftheCFDApproach
giveninTableII.Lactosepowderwassimulatedinthisstudy
in order to compare the computational results with
experimentalresultsfrom(41)forvalidationpurposes.While The CFD method developed and reported in this study has
thereallactosepowdercomprisesofparticleswithadistribu- beenvalidatedbycomparisonwithentrainmentexperiments
tion of particle sizes, in this study a mono-dispersed powder performed by Tuley et al. (41). In Tuley’s study, a constant
was simulated. Lactose particles are usually used as diluent, pressuregradientdpwasappliedattheoutletofthreedifferent
dt
alsoknownascarrierparticles(42,43)andhavehighaverage entrainment geometries. Each geometry permitted entrain-
diameter compared to the drug, excluding the fines. In this ment by a different combination of shear and gas assisted
studyde-agglomerationwasnotsimulated.Althoughasimpli- mechanisms. The entrainment geometries were transparent
ficationforthepurposesofthecurrentstudy,itwasassumed and a video camera was used to record the entrainment of
thatthedrugparticleswouldde-agglomerate fromtheir lac- particles.InordertovalidatetheCFDapproachchoseninthis
tose carrier during or after entrainment. Delivery of active study, one of these entrainment geometries was modelled
drug is therefore assumed to be proportional to the amount usingtheCFDapproachdescribedinSections2.3.1–2.3.3.
ofentrainedlactosepowder. Forthegeometryselected(see Fig.5) theinitial entrainment
Kopsch,MurnaneandSymons
meshIndependenceStudy
TheaccuracyandthecomputationalcostsofCFDsimulations
areinfluencedbythetypeofchosenmesh(44).Ingeneral,a
high mesh density results in more detailed and accurate re-
sults, but also higher computation time. It is important to
verifythatthenumericalresultsareindependentofmeshres-
olution.Smallchangesinthemeshshouldnotresultinlarge
changes in the results. A mesh-independence study may be
conducted by plotting one variable for different meshes.
Fig.6 shows the result ofa m(cid:4)esh-independenc(cid:5)e study:For a
particular geometry a*;b*;c*;d*;e* ¼ð0:5;0:5;
0:5;0:5;0:5Þ,thenumberofcellsusedtomeshthegeometry
wasvaried.Itappearsthatmesheswithmorethan3750cells
aresufficienttoachievemeshindependence.
ExplorationofDesignSpaceforDPIOptimization
Beforethedesignspacecouldbeexploredtwoparameters
had to be set. First, the relative weightings w and w of
A B
the cost functions C and C had to be set. Second, the
A B
value of x in the ideal drug release profile M (x) had to
0 3
be set (see Eq. 7). Initially, the weightings were set to
w=(w ,w )=(1,1),i.e.bothcostfunctionswereweighted
Fig.5 Comparisonofexperimentaldata(Tuleyetal.(41),Fig.6)withnu- A B
equally. However, C tended to be (much) smaller than
mericalsimulation.Theinletistherightleg,theoutletistheleftleg.Tuleyetal. A
C . This implies that for the modelled (pierced blister)
(41)filledlactosepowder(16%fines)intotheentrainmentdevice. B
geometry, objective A (similar drug release for different
patients) is comparatively easier to achieve than objective
B (similarity to a desired drug release profile).
mechanismisgas-assisted,butittransitionstosheardrivenas Consequently, the weighting was changed to w=(40, 1)
evacuationprogresses.
in a second optimization process. These weighting values
Figure5showsimagesofanentrainmentdevicefordifferent
imply that most emphasis is put on achieving objective A.
timesaslactosepowder(16%fines)wasentrainedwhenapres-
In this case objective B is (almost) ignored.
suregradientofdp=−30kPas−1wasappliedattheoutletofthe In addition, two different values for x , the duration of
dt 0
device (41). For comparison, the numerical results (EE CFD) drug release in the ideal release function M (x), were ex-
3
are shown next to the experimental results. Tuley et al. (41) amined(seeEq.7).Avalueofx =0.2givesanidealearly
0
recordedthefirstframeatt ¼0:033s,andthereforeallCFD bolus drug release, while a higher value of x =0.5 gives a
0
modellingtimepointsaregivenrelativetothistime.Theplotsof morecontinuousreleaseprofile.Asanexample,theresults
streamlines shown in Fig. 5 have been compiled from the air
velocitysolution(EECFD).Thestreamlinesarecurvesthatare
tangent to the velocity vector of the air at any given time.
Typically, streamlines converge when the speed of the fluid
increases.Therefore,plotsofstreamlinesmaybeusedtoana-
lysehowtheflowisaffectedbythepresenceofthedrugpowder:
Initially, the volume fraction of drug α is 0.49, see Table II.
Oncedrugismovedawayfromthepowderbedsurface,athin
channel of air is shaped. The streamlines converge towards
eachotherinthischannel.Thus,airflowspreferablythrough
this channel. The results from the simulation appear to be
similartotheexperimentalresultswiththeformulatedpowder
beingreleasedoversimilartimescalesandinasimilarprocessin
themodelasintheexperimentalobservations. Fig.6 Meshindependencestudy.
OptimizingdryPowderInhalers:Methods+Results
Fig. 7 Release profile and costs for different values of the length of the compartment end e*. All other parameters were kept constant:
(a*,b*,c*,d*)=(0.76,0.26,0.20,0.38).(a)Examplereleaseprofilesforfourdifferentvaluesofe*.Theinhalationprofilesappliedwerethoseofaweaker
patient(Q2(t)).(b)ValueofthecostfunctionsCAandCBase*isvaried.
for the scaled release profile and the corresponding cost OptimizationResults
functionsarepresentedinFig.7whenjustasingleparam-
∗
eter, the length of the compartment end e , was varied, Optimization of the entrainment geometry of a DPI was
while all other parameters were kept constant. Fig. 7a attempted to achieve two objectives. One objective was to
shows that early delivery of drug is achieved for small deliver drug to similar parts of the airways for different pa-
∗
values of e , while a more continuous delivery is achieved tients(objectiveA).Thesecondobjectivewastodeliverdrug
∗
forlargervaluesofe .Fig.7b indicatesatrade-offbetween toaspecificlocationintheairways(objectiveB).Optimization
∗
costs C and C : A high value of e results in a continuous was conducted for four different combinations of the
A B
drug delivery and hence a small value of C (for x = 0.5). weightings w=(w , w ), of the cost functions (C and C ),
B 0 A B A B
Howeverinthisregime,oflongerdurationdelivery,differ- andthedeliveryduration,x (seeFigs.8,9and10).Ingeneral,
0
ences between the emission profiles for patients 1 and 2 theoptimizedgeometrydependsontheweightingsofthecost
become more significant and so C increases. functionwandthedurationx .Figs.8,9and10showopti-
A 0
The results show that powder entrainment can be influ- mizationresultsfordifferentvaluesofwandx .
0
encedbychangingdesignparameters.Inparticular,abetter Figure8showstheoptimizedgeometry(a)andthecorre-
early bolus delivery may be achieved. This phenomenon is sponding scaled drug release profiles (b) when the cost func-
consistentwithexperimentalfindings.Forexample,Tuleyetal. tionswereequallyweightedw=(1,1),i.e.C =C +C ,and
tot A B
(41) compared different entrainment devices and found that anearlybolusdelivery(x =0.2)wasassumed.The‘highflow
0
boththefluidizationprocessandthetimescaleoffluidization rate’profilecorrespondstothedrugreleaseprofileforpatient
isinfluencedbythechoiceofgeometry. 1 with average inhalation flow, Q (t). The ‘low flow rate’
1
Fig.8 C=CA+CB,x0=0.2,the
optimalsetofdesignvariablesis
x=(0.68,0.10,0.64,0.47,0.08).
(a)Theoptimalgeometryandthe
volumefractionαattimet=0.2s.
(b)Scaleddrugreleaseprofilefor
highandlowflowrate.
Kopsch,MurnaneandSymons
Fig.9 C=CA+CB,x0=0.5,the
optimalsetofdesignvariablesis
x=(0.55,0.0,0.79,0.34,0.93).
(a)Theoptimalgeometryandthe
volumefractionαattimet=0.2s.
(b)Scaleddrugreleaseprofilefor
highandlowflowrate.
profile corresponds to the drug release profile for patient 2, toachievemorecontinuousdrugdeliveryhighervariationsin
withlowerinhalationflowQ ðtÞ.Similarly,Fig.9showsthe thedeliveryfordifferentpatientshavetobeacceptedforthis
2
optimized geometry (a) and the corresponding scaled drug typeofgeometry(asshowninFig.9b).
releaseprofiles(b)forthesameweighting,butamorecontin- Fortheexamplegeometrychoseninthisstudy,nodevice
uousdrugdelivery(x ¼0:5).Comparingtheoptimizedge- was found that effectively delayed the drug delivery in the
0
ometries for these two cases suggests that a continuous drug scaledrepresentation.Therefore,onlyabolus-deliveryatthe
deliveryisachievedbyincreasingthelengthofthecompart- beginning of the inhalation cycle or a continuous delivery
mentende.Drugthatisplacedwithinthecompartmentend from the beginning of the inhalation cycle are possible.
isslowlyentrainedbyashearingactionduringinhalation. Developing a passive DPI geometry that allows a delay in
WhentheC andC weightingsofw=(40,1)werechosen timing of the drug delivery will be the subject of future
A B
(i.e. emphasis was put on objective A to achieve maximum research.
similarityofdrugreleaseinthesameinhaledvolumeelements
for patients with different inhalation manoeuvres), approxi-
mately the same optimized geometry was discovered, inde- CONCLUSIONS
pendent of the timing value x that was chosen. Figure 10
0
shows theoptimizedgeometry andthecorrespondingscaled Inthiswork,acomputationalmethodhasbeendevelopedand
drug release profiles when weightings of w = (40,1) and a applied for in silico optimization of a simple DPI geometry.
timing valueof x = 0.2 was chosen.However, itwas found Particularly,twocostfunctions thatquantifyperformance of
0
thatverysimilarresultswerediscoveredforatimingvalueof DPIs and a dynamic CFD boundary condition were devel-
x = 0.5 and thus only one figure ofthis typeisshown.The oped. WehaveshownthatCFDdesignoptimizationcanbe
0
scaled drug release profiles indicate that good agreement of appliedtoDPIdesign.Validationisanimportantissuewhen
releaseprofilesofdifferentpatientswasachievable,butopti- consideringnumericalsimulations.Duetothelackofexperi-
mizationobjectiveB,i.e.achievingatargetreleaseprofile,was mentalresults,validationwasonlypossibletoalimitedextent
effectivelyignoredduetothelowrelativeweighting.Inorder in this study. Nevertheless, the computational approach
Fig.10 C=40CA+
CBx0=0.2andx0=0.5,the
optimalsetofdesignvariablesis
x=(0.76,0.29,0.21,0.38,0.0).
(a)Theoptimalgeometryandthe
volumefractionαattimet=0.2s.
(b)Scaleddrugreleaseprofilefor
highandlowflowrates.