Stochasticity in Primary Nucleation: Measuring and ... - ACS Publications

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Stochasticity in primary nucleation: measuring and modelling detection times Giovanni Maria Maggioni, and Marco Mazzotti Cryst. Growth Des., Just Accepted Manuscript • Publication Date (Web): 31 May 2017 Downloaded from http://pubs.acs.org on June 6, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Crystal Growth & Design is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

Stochasticity in primary nucleation: measuring and modelling detection times Giovanni Maria Maggioni and Marco Mazzotti∗ Separation Processes Laboratory, ETH Zurich, Zurich E-mail: [email protected] Phone: +41 44 632 24 56. Fax: +41 44 632 11 41 Abstract

1

2

30th May 2017

3

In this work, stochastic nucleation has been investigated analysing a large set of

4

literature data on the isothermal crystallisation of p-Aminobenzoic acid in acetonitrile,

5

isopropanol, and ethyl acetate in 1.5/1.8 mL batch reactors. We investigated the

6

basic statistical hypotheses used to model the stochastic process of nucleation and

7

the criteria necessary for their fulfilment: based on a sound mathematical framework,

8

we have developed and implemented a statistical analysis which allows to assess the

9

statistical representativity of experimental data and their inherent uncertainty.

10

A physical model of the process has been developed including primary nucleation

11

and growth of crystals until detection. Based on the statistical analysis developed,

12

the model parameters have been estimated from the measurements, accounting for

13

stochastic nucleation and conditions for detection, according to two alternative meth-

14

ods.

15

In light of such results, this work explores two further aspects in details: the choice

16

of the appropriate stochastic process to describe nucleation, and the sensitivity of the

17

measurements to operating parameters, with a special emphasis on supersaturation

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

18

and temperature. Finally, we give recommendations how to measure and interpret

19

stochastic nucleation data.

20

1

Introduction

21

Nucleation is the formation of entities of a new phase from a parent phase 1,2 ; in the case of

22

crystallisation, new particles of a species dissolved in a solvent are formed from its homo-

23

geneous solution upon creation of supersaturation, i.e. upon reaching a solute concentration

24

above its thermodynamic limit, i.e. its solubility. To understand and describe nucleation

25

we need to characterise its thermodynamics and its kinetics. Particularly kinetics, i.e. more

26

specifically an expression for the rate of nucleation as a function of temperature and solute

27

concentration, is essential to be able not only to describe nucleation in the context of a

28

population balance model of the whole crystallisation process, but also to understand its

29

underlying physics.

30

Measuring nucleation kinetics has always been and still is a challenge, first because new

31

crystals can be detected only after they have grown to a certain minimum size (dependent on

32

the detection technique) and because nucleation experiments seem to exhibit an intrinsic lack

33

of good reproducibility. The former issue makes it necessary to study nucleation coupled with

34

crystal growth, unless one can make the assumption of fast growth and negligible growth time

35

between nucleation and detection. The latter issue has been tackled by using smaller and

36

smaller devices (down to the use of microfluidics) where conditions can be better controlled,

37

and by running a large number of experiments so as to build an adequate statistics. Such

38

approach has not been entirely successful, in the sense that on the one hand nucleation times

39

have been spread over even a larger interval than in larger crystallisers 3,4 , and on the other

40

hand the possibility of measuring nucleation rates at a small scale and then extrapolating

41

such information to larger crystallisers appears still problematic 5,6 .

42

Today, it is accepted that the formation of new nuclei is a stochastic phenomenon 3,7,8 ,


Page 2 of 38

Page 3 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


43

governed by for instance the law of the occurrence of rare events, as a consequence of the

44

existence of an energy barrier in the formation of new nuclei 1,2,9,10 and of the fact that the

45

first nucleus or the first nuclei in a suspension are by definition individual events, to which a

46

treatment based on averaging their properties cannot be applied (as in the case of chemical

47

reactions).

48

The stochastic nature of nucleation has long been recognised 11 and in the last years it

49

has been studied with increased motivation and efforts both theoretically 3,4,8,12–14 and exper-

50

imentally 4,5,7,15–17 . More powerful computational tools and in-situ experimental instruments

51

have triggered such development, together with the availability of medium-throughput tools,

52

e.g. the Crystal16 (manufactured by Technobis Crystallization Systems), that make it pos-

53

sible to carry out a large number of repetitions of the same identical experiment.

54

The abundance of new experimental data and the stochastic approach to the study of nuc-

55

leation have helped in formulating new theoretical interpretations and models 3,4,6,8,12,18–22 ,

56

but, at the same time, have yielded some new and important challenges for both experi-

57

mentalists and modellers 17,23–28 . First of all, each set of experiments analysed should be

58

truly representative of the stochastic process. On the one hand, it should be possible to

59

reproduce experiments in a reliable way, so that the statistics extracted from the data effect-

60

ively belongs to the same underlying distribution. On the other hand, enough measurements

61

must be collected to properly reconstruct the actual probability distribution function from

62

the experiments, as it is the case for every statistical phenomenon. Were the data not rep-

63

resentative of the parent statistics, their analysis, and the results of such analysis, would

64

be meaningless. A second issue, strictly correlated to experiments reproducibility, concerns

65

the sensitivity of the system to small variations of initial and boundary conditions (always

66

occurring within different repetitions of a theoretically identical set of experiments). A sys-

67

tem should not be sensitive to such variations in a strong way, otherwise it would not be

68

possible to discriminate between variations due to the stochastic behaviour and those due

69

to experimental uncertainties. The latter, even if fundamentally different in nature, can



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

70

produce results that appear to be stochastic. Despite their similarity, though, these types of

71

variabilities have a different physical nature and interpretation.

72

A third issue is that of how to estimate nucleation kinetics from detection time meas-

73

urements characterised by an inherent stochasticity and affected by the sensitivity of the

74

detection instrument. The literature shows how additional assumptions about detection and

75

growth have to be made 4,5,29,30 , how nucleation and growth kinetics have to be estimated at

76

the same time 6,29 , and how the stochasticity and the variability of nucleation rate measure-

77

ments propagate to the model parameters belonging to the nucleation rate expression 6,28 .

78

The aim of this work is to study the issues above by analysing one of the most extensive

79

set of data available to date, namely those referring to the nucleation of α p-Aminobenzoic

80

Acid (p-ABA) in three different solvents, i.e. acetonitrile (ACN), 2-propanol (IPA), and

81

ethyl acetate (EtOAc), which were recently published by Sullivan et al. 17,31 . Our primary

82

goal is that of analysing the statistical features of this wealth of data and of understanding

83

the limitations of using a finite set of data to extract reliable, supersaturation-dependent

84

nucleation kinetics (one set of rate parameters for each solvent of course; see Sections 3 and

85

4 of this work). In doing so, one needs to estimate also the parameters in a growth rate

86

expression, since crystals are observed not when they form, but after a certain extent of

87

growth, as discussed elsewhere 6 . After providing the theoretical statistical background on

88

which our statistical analysis is based (Section 2), we present a physical model aimed at

89

describing the physical process of nucleation and crystallisation from solution (Section 3).

90

Then, we report and discuss the results of applying our model to the experiments of Sullivan

91

et al. 17 (Section 4). In order to consider all possible reasons for discrepancies between

92

the model and the data, we comment on the choice of the stochastic process and on the

93

sensitivity of the system to operating parameters (Section 5), before summarising the results

94

and drawing some conclusions (Section 6).


Page 4 of 38

Page 5 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


95

2

Statistical analysis of nucleation experiments

96

2.1

97

Detection time experiments have been repeatedly observed and reported to be the result of a

98

stochastic process, as illustrated in Figure 3 in the work of Sullivan et al. 17 and in Figure 2 in

99

our earlier paper 6 , that uses data from another paper 5 . A proper analysis and quantitative

Fundamental hypotheses and statistical preliminaries

100

interpretation of such experiments rest on three fundamental assumptions:

101

H.1) the experiments are carried out under identical and reproducible experimental conditions;

102

103

H.2) the detection times are stochastic events, which sample one and the same cumulative

104

probability function, which depends on the chemico-physical properties of the system;

105

H.3) the experiments are in a sufficient number to form a representative sample of the parent statistics.

106

107

Some basic definitions and concepts are needed, to provide a sound understanding of these

108

hypotheses and of their implications. Let us consider a set of Mo observed, independent

109

detection time measurements tD , extracted from M ≥ Mo independent experiments, such

110

that: D D T D D tD = tD : tD 1 , t2 , ..., tMo 1 ≤ t2 ≤ ... ≤ tMo

111

(1)

The measured empirical Cumulative Distribution Function (eCDF), F, associated to tD is: F = [F1 , F2 , ..., FMo ]T : Fj =

j , j = 1, 2, ..., Mo M

(2)

112

Mo is the number of experiments where crystals were detected, while (M −Mo ) is the number

113

of experiments where no detection occurred. In the following, tD , tN , F , without index, will be



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 38

114

used to denote generic quantities, whereas the indexed quantities will correspond to specific

115

experiments. The statistical nature of the detection times is mathematically expressed by

116

saying that tD is sampled from a distribution, whose cumulative distribution function is P (t),

117

in symbols:

tD ∼ P

(3)

118

This expression should be read as “tD is distributed according to the cumulative distribution

119

P ”. Let us also define:

P = [P1 , P2 , ..., PMo ]T : Pj = P (tD j ), j = 1, 2, ..., Mo

(4)

120

The figures mentioned above, namely Figure 3 in Sullivan et al. 17 and Figure 2 in Maggioni

121

and Mazzotti 6 , plot the empirical distribution F vs the detection time vector tD . In case of

122

a model-based description of the data also the cumulative distribution function P is plotted

123

vs the vector tD . F is an approximation of P, and the distance between F and P is measured

124

by means of the Chebyschev norm, D, also used in the Kolmogorov-Smirnov test:

125

D = kF − Pk∞ = sup |Fj − Pj |

(5)

1≤j≤Mo

126

The Glivenko-Cantelli theorem 32 proves the uniform convergence in probability of F to P: a.s.

lim D = 0

Mo →∞

(6)

127

Note that the Glivenko-Cantelli theorem in its general mathematical statement uses the limit

128

M → ∞, since in theory one can assume that the experiments run as long as needed in order

129

to let the system produce an observation. Here, though, to be consistent with this physical

130

limitation, we have chosen to use Mo , since M approaches infinity if Mo does, while the


Page 7 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


131

opposite is not true. Eq. (6) assures that, given a sufficient number of successful experiments,

132

i.e. experiments where detection was observed, the probability built from experiments using

133

Eqs. (1) and (2) represents a reliable approximation of the real probability, whatever this

134

is, which means that the sequence D approaches 0 as Mo → ∞ almost surely, i.e. with

135

probability 1. Note that D depends not only on the number of experiments, M , but also

136

on the number of observations, Mo , and that when Mo → ∞ also M → ∞. The rate of

137

convergence of this sequence will be discussed in Section 2.2, in relation to the fulfilment of

138

H.3.

139

2.2

140

Since we aim at modelling a stochastic behaviour based on experimental evidence, it is essen-

141

tial to investigate the validity of hypothesis H.3, i.e. if and when the number of experiments

142

M is large enough to constitute a statistically representative sample of the parent population.

143

The Glivenko-Cantelli theorem, Eq. (6), guarantees uniform convergence in probability of

144

the experimental cumulative distribution F to the underlying distribution P, but does not

145

provide any estimation of the rate of convergence, hence it cannot be used to test H.3. A

146

minimum rate of convergence of Eq. (6) is estimated by the Dvoretzky-Kiefer-Wolfowitz

147

inequality 33,34 (DKW), which places probabilistic bounds (DKW-bounds) on the maximum

148

value of D, i.e. the distance between F and P, defined in Eq. (5), and reads:

Representativity

P (D ≤ ε) ≥ 1 − 2e−2M ε

2

(7)

149

where ε (> 0) represents the upper bound of D, M is the number of experiments, and

150

P (< 1) is the probability that D is smaller than ε. From Eq. (7), a lower bound on the

151

number of points M to achieve a given accuracy ε with a given probability P is estimated



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

152

Page 8 of 38

as: 1 M ≥ 2 ln 2ε

2 1−P

(8)

153

and a lower bound on the maximum value of D for a given number of points M with a given

154

probability P as: s ε≥

1 ln 2M

2 1−P

(9)

155

Eq. (8) is useful for the design of experiments, namely to choose the number of measurements,

156

M . Eq. (7) and (9) are useful to analyse the quality of a given set of M measurements.

157

Quite significantly, these results are independent of the underlying statistics, i.e. they do

158

apply to any stochastic process, e.g. to the Poisson processes. One should also note the

159

strong non-linear dependence of M on P and ε, which leads to a very rapid increment of M

160

if high accuracy (small ε) and high reliability (large P) of F are required.

161

3

162

Crystal nuclei cannot in general be observed directly when they form, since they are too

163

small for most measurement techniques. They can be measured only after growing to a de-

164

tectable, finite size. For this reason, one speaks correctly of measurements of detection time,

165

which is defined as the time elapsed from the establishment of supersaturation to when some

166

measurable property of the system reaches and/or overcomes a threshold value. The quant-

167

ities which can be monitored on-line are usually the turbidity (turbidimetry/nephelometry),

168

the infra-red absorbance (ATR-FTIR), or the reflectance (FBRM) of light in the suspension.

169

The estimation of nucleation kinetics is performed by correlating these measured quantities

170

to properties of the crystal distribution (typically moments of the distribution, see below

171

Eq. (19)) and/or of the solution, for example:

Physical and mathematical model


Page 9 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


172

• the number of crystals (estimated from FBRM measurements);

173

• the crystals’ average size (FBRM);

174

• the crystals’ volume fraction (from turbidimetry and FBRM measurements);

175

• the total crystal surface per unit volume of the suspension (FBRM, turbidimetry)

176

• the solute concentration (measured by ATR-FTIR).

177

In order to interpret the detection time measurements, a model that describes how the

178

system evolves from nucleation to detection is necessary. Such a model requires additional

179

assumptions, namely:

180

H.4) the detection time is formed by the additive contribution of nucleation and growth time.

181

The nucleation time is the time elapsing between the establishment of supersaturation

182

and the formation of the first nucleus, whereas the growth time is the time elapsing

183

between the nucleation time and the detection of crystals;

184

185

186

187

H.5) the observed stochastic behaviour of detection time is entirely due to primary nucleation, whereas crystal growth is treated as a purely deterministic phenomenon; H.6) the first nucleus in the system is the only stochastic event, and follows a well-defined stochastic process.

188

Mathematically, one can envisage the nucleation time tN , the growth time tG , and the

189

detection time tD as the (model) output of the operators H N , H G , and H D , respectively.

190

These operators determine the evolution of the crystallising system having as input the

191

initial concentration, c0 , the system size w (which can be the mass of solvent, or the volume

192

of the solution, and is assumed in this work to be constant for simplicity, but without loss

193

of generality), the temperature evolution, T (t), and the detection condition, α, hence they



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

194

Page 10 of 38

return the relevant characteristic times:

tD = H D [c0 , T (t), w, α],

(10)

tN = H N [c0 , T (t), w],

(11)

tG = H G [c0 , T (t), w, α]

(12)

196

tD = tN + tG

(13)

197

It is thus necessary to find appropriate functional forms for H N and H G , whereas H D can be

198

obtained from their simple summation. Our goal has been to create a physically consistent

199

and meaningful, though complex, model which not only minimises the number of parameters

200

required at its implementation, but also assigns them physical meaning.

195

with:

201

In agreement with H.4 and H.5, after the formation of the first nucleus, the system

202

evolves in a deterministic way, as described by a 1D population balance equation (PBE) for

203

the evolution of the crystal size distribution (CSD), f (t, L) 35,36 . As to the CSD definition,

204

wf dL is the number of particles in the suspension with size between L and L+dL; moreover,

205

the moment of order i of f is defined as: Z

206

φi (t) =

∞

Li f (t, L)dL

(14)

0

207

The evolution of the properties of the CSD is described also by the Equations of Moments

208

(EoM), that can be derived from the PBE. For a standard batch crystalliser, whose temper-

209

ature evolution is followed exactly thanks to an ideal thermostat, in case of size-independent


Page 11 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

210


growth, the EoM read: dφi = Jδi0 + iGφi−1 dt φ0 (tN ) =

i = 0, 1, 2, ... 1 w

φi (tN ) = 0

(15) (16)

i>0

(17)

211

where G is the growth rate, J is the nucleation rate, and δi0 is 1, if i = 0, and 0 otherwise.

212

Only a finite set of these equations needs solving, namely always for i = 0, 1, 2, 3, possibly

213

also for i = 4 for certain detection conditions. The set of equations (15) is coupled to the

214

mass balance, that provides the solute concentration c, and reads:

215

c(t) = c0 − kv ρc φ3 (t)

216

The initial conditions of the model equations are assigned at t = tN , hence eqs. (14), (15)

217

and (18) are defined for t ≥ tN . Note that according to eqs. (16) and (17) it is assumed that

218

only one crystal of negligible size is present after the first nucleation event at the beginning

219

of the deterministic process.

(18)

220

The integration of the EoM is performed until a stop condition is fulfilled, and the time

221

associated to it is the model detection time tD mod . Such condition is typically given in terms

222

of a function of one or more moments of the distribution reaching a specific value, i.e.

D D D Φ(φ0 (tD mod ), φ1 (tmod ), φ2 (tmod ), φ3 (tmod ), ...) = α

(19)

223

where the function Φ depends on the physical principle behind the measurement and α is

224

system- and instrument-specific threshold value.

225

Let us now consider the constitutive equations that describe specific properties of the

226

system, e.g. the rates of crystal nucleation and growth. First of all, the supersaturation

227

S is defined as the ratio between the actual solute concentration, c, and the temperature 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

228

dependent solubility, c∗ (T ): S=

229

230

Page 12 of 38

c c∗ (T )

(20)

The solubility c∗ (T ) can be expressed in terms of T for instance through:

231

q 1 c∗ (T ) = q0 exp − T

232

with q0 , q1 being two constants; here c∗ and q0 are in gram of solute per gram of solvent

233

[g/g], while T and q1 in [K]. The values of the solubility parameters, (q0 ; q1 ), have been

234

obtained from fitting Eq. (21) to the data provided by Svard et al. 37 ; they are (328; 2, 510),

235

(224; 2, 380), and (3.1; 1, 078), for ACN, IPA, and EtOAc, respectively. Here it is worth

236

making a comment, namely that our analysis focuses on nucleation triggered by temper-

237

ature changes (cooling crystallisation). However, a similar analysis could be repeated for

238

evaporative crystallisation (where w would be a function of time ) and for antisolvent crys-

239

tallisation (where w = w(t) and c∗ would be function of T and of the antisolvent fraction).

240

241

(21)

Through S, we define the growth rate G, which appears in Eq. (14), according to the birth-and-spread growth model:

242

K1 K2 2/3 1/6 G(S, T ) = K0 exp − (S − 1) (ln S) exp − 2 T T ln S

243

This choice is justified by the conclusion drawn in a recent paper 38 and by the fact that the

244

range of supersaturation values explored by Sullivan et al. 17 was limited to between 1.05 and

245

1.30, thus not justifying considering multiple growth mechanisms.

246

247

248

(22)

We assume that the stochastic process postulated by H.6 for nucleation is a Poisson process 3–5,8,11 , hence P (t) has the following exponential form:

P (t) = 1 − exp (−Λ(t))


(23)

Page 13 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


249

where Λ(t) represents the expected number of primary nuclei formed by time t, and is in

250

general given by the following integral function: t

Z

wJp (s)ds

Λ(t) =

251

(24)

0

252

where Jp (s) is the primary nucleation rate prevailing at time s. The lower limit of the

253

integral is zero, since we set as initial time of the process the moment when the chosen

254

nominal supersaturation (S > 1) is attained in the system, and because the nucleation rate

255

is identically zero for saturated and undersaturated conditions. It is worth noting that in

256

Eq. (24) what determines the value of Λ(t) is the product of Jp and w, rather than their

257

individual values. The instantaneous primary nucleation rate Jp depends on time through

258

the corresponding concentration and temperature values. In agreement with the conclusions

259

of Sullivan et al. 17 , we assume that it follows the Classical Nucleation Theory 1 (CNT) given

260

by:

261

A1 Jp (S, T ) = A0 exp − T

B S exp − 3 2 T ln S

(25)

262

where A0 , A1 and B are constant. Other mechanisms, e.g. the two-step mechanism proposed

263

by Vekilov 39 , might be included here if there were evidence that they were dominant.

264

While in general the nucleation rate J in Eq. (15) could be the sum of primary and

265

secondary nucleation, it is customary to consider only primary nucleation, i.e. J = Jp to

266

describe detection time experiments. We do the same in this paper, thus being consistent

267

with previous literature 5,6,17,28,30 . In this case, we consider as detection condition a threshold

268

value on the average number weighted crystal size, namely:

269

Φ=

φ1 = αL φ0

(26)

270

with αL = 10 µm in the computations below. Φ is representative of a (number weighted)

271

average crystal size, and the chosen value of αL is consistent with other papers in the liter-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 38

272

ature 28 . The set of model parameters consists then of ω = [A0 , A1 , B, K0 , K1 , K2 ]T , and has

273

to be assigned to solve the model equations.

274

4

275

In the previous section, we have established the mathematical model used to describe the

276

experimental system, which consists of equations (15)-(18), with the detection condition

277

given by Eq. (19). The operating conditions, i.e. temperature T (t), system size w, and initial

278

solute concentration c0 , are the input parameters. If the system is operated at isothermal

279

conditions, Eq. (24) simplifies to Λ(t) = wJp t, K0 and K1 in Eq. (22) can be lumped into

280

a single parameter, i.e. K = K0 exp(−K1 /T ), and A0 and A1 in Eq. (25) can be combined

281

into A = A0 exp(−A1 /T ). In this section, we discuss how ω can be estimated by fitting

282

the experimental results. It should be noted that we express the system size w as mass of

283

solvent, and thus the units of Jp are [# kg-1 s-1 ].

Results

284

The experimental results considered in this work have been reported by Sullivan et al. 17 ,

285

who carried out isothermal crystallisation of p-Aminobenzoic acid in vials with a volume

286

of 1.5 or 1.8 mL, at a stirring rate of 700 rpm, in three different solvents (Acentonitrile,

287

ACN, 6 series of experiments at different supersaturation levels; Isopropanol, IPA, 8 series;

288

Ethyl-acetate, EtOAc, 7 series). The necessary input values consist, for each experimental

289

series in a specific solvent, of the set of measured detection times, i.e. the vector tD of Eq.

290

(1), and the corresponding eCDF, i.e. the vector F of Eq. (2). Since tN is equal to tD minus

291

a deterministic constant (see eq. (13) and hypothesis H.5), tN is distributed as tD , i.e.

292

tN ∼ P

293

hence the eCDF is representative of the distribution not only of detection times, but also of

294

nucleation times.

295

(27)

There are two different possible metrics for the distance between the experimental data 14 ACS Paragon Plus Environment

Page 15 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


296

and the simulation results, i.e. the error. Since the available experimental data are the

297

vectors tD and the eCDF F, such errors can be calculated based on either the former or the

298

latter. The first uses for each series k the distance ejk between each detection time forecast

299

by the model and its experimental counterpart; such distance is defined as the difference

300

between experimental and model values, relative to their geometric mean:

301

D tD jk − tmod,jk ejk = q D tD jk tmod,jk

302

The use of the weight, namely the reciprocal of the geometric mean in this case, guarantees

303

that the fitting procedure is not biased in favour of longer detection times. The total error

304

for each individual series k is thus computed as the euclidean norm Ek of the vector ek ,

305

consisting of the ejk associated to the k th series:

306

Ek = kek k2 =

sX

(28)

e2jk

(29)

j

307

Finally, the total error for each solvent, Es , is given by the euclidean norm of the vector E,

308

constituted of the errors Ek associated to the different series of experiments, divided by the

309

number of series Ns for that solvent:

310

1 1 Es = kEk2 = Ns Ns

sX k

Ek2

1 = Ns

sX X k

e2jk

(30)

j

311

The second metric is based on cumulative probabilities instead of detection times, and

312

uses the Chebyschev norm (supremum norm) introduced in Eq. (5). For each series k

313

belonging to a specific solvent we define Dk as:

314

Dk = kFk − Pk k∞

315

where Fk and Pk are the distribution measured or calculated at the conditions of series k. 15 ACS Paragon Plus Environment

(31)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

316

The Dk values substitute the Ek values of the first metric; how to obtain an overall index of

317

the model accuracy for the specific solvent, i.e. the p-value ps , will be discussed in Section

318

4.3. Model parameters are estimated by optimising either Es or ps . The two approaches lead

319

in principle to different parameters and we aim at understanding which method is best and

320

why. We present and discuss the methods, and their results in the following in detail. All

321

c optimisations have been performed using the Matlab routine fminsearch, a derivative-free

322

method.

323

Before proceeding with the discussion of the parameter estimation procedure, it is worth

324

making three points. First, p-ABA exhibits rather fast growth kinetics, which implies very

325

short growth times and detection occurring very soon after nucleation. Under these con-

326

ditions and for isothermal experiments (not linear cooling experiments, as analysed in our

327

previous contribution 6 ) the nucleation rate and the growth rate parameters are not correl-

328

ated, as we could demonstrate by estimating the nucleation rate parameters repeatedly for

329

different assigned values of the growth rate parameters; by doing so, we have always obtained

330

very similar nucleation rate parameters. Second, it is worth keeping in mind that the nuc-

331

leation rate parameters are strongly dependent on the threshold value used in the detection

332

condition of Eq. (26). We have chosen the most sensible value for such threshold, namely a

333

value in accordance with earlier papers 17,28 . Finally, we have estimated confidence intervals

334

for the parameters according to a heuristic procedure, i.e. we have repeated the estimation

335

of the parameters from randomised initial guesses and looked for the minimum of the errors,

336

using a genetic algorithm. The minimum found by the optimiser, i.e. the genetic algorithm,

337

was assumed to be also the global minimum; the parameters associated to the local minima

338

differing not more than 1% from the global minimum were used to estimate the standard

339

deviation, which was used as confidence interval for the parameters.


Page 16 of 38

Page 17 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


340

4.1

Nonlinear least-squares method (NLSM)

341

Minimising Es defined by Eq. (30) implies solving a non-linear least squares problem, which

342

entails establishing an iterative procedure, at each n-th step of which a set of paramet-

343

ers, ω (n) , is computed as current guess. With these, tD mod,jk in Eq. (28) can be calcu-

344

lated as follows. First, the experimental Fjk value associated to the measured tD jk value is

345

N used to compute tN mod,jk from the isothermal version of Eq. (23), thus obtaining tmod,jk =

346

−(ln(1 − Fjk )/(wJk )), where Jk is the nucleation rate at the conditions (supersaturation and

347

temperature) of the series k. Then, Eqs. (15)-(19) are integrated from tN mod,jk to when the

348

detection condition (Eq. (19)) is fulfilled: this establishes the time corresponding to the

349

requested value tD mod,jk .

350

The experimental results have been described with the model above. All experiments in

351

the same solvent have been fitted together thus obtaining the three sets of parameters repor-

352

ted in Table 1. The corresponding errors for each series of experiments at the corresponding

353

value of the supersaturation are reported for all twenty-one series in Table 2. Table 1: Nonlinear least-squares method. Estimated values of the model parameters for the three solvents. Model parameters, ω A [# kg−1 s−1 ] B × 10−5 [K3 ] K × 107 [m s-1 ] K2 × 10−3 [K2 ]

ACN IPA EtOAc 0.94 ± 0.1 0.44 ± 0.1 0.42 ± 0.05 2.81 ± 0.1 4.45 ± 0.2 4.11 ± 0.2 1.78 ± 0.5 9.67 ± 2 3.37 ± 0.7 8.61 ± 0.2 1.58 ± 0.1 0.92 ± 0.2

354

355

The comparison between experimental data and estimated results is illustrated in Figure

356

1, where the cumulative distribution functions obtained in each series of experiments, i.e.

357

the points (tD , F), are compared with the curves calculated using the model, i.e. (tD mod , P).

358

In each figure, the supersaturation level increases from lower right to upper left, i.e. going

359

from longer detection times to shorter ones. 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Comparison between experimental measurements and simulation results calculated using the parameters estimated through the nonlinear least-squares method. The experimental points of Sullivan et al. 17 are represented by grey diamonds, while the calculated distributions are plotted as continuous lines, from black to red, indicating increasing supersaturation. The abscissa is in logarithmic scale. 360

Let us consider the experimental points first. It is rather apparent that some of the

361

experimental curves exhibit a rather regular behaviour, e.g. those in ACN in general and a

362

few in IPA and EtOAc. On the contrary, a few curves have a shape rather different from

363

the exponential Poisson distribution, e.g. in IPA at the highest supersaturation. In other 18 ACS Paragon Plus Environment

Page 18 of 38

Page 19 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


364

cases the experimental behaviour appears to be counter-intuitive, e.g. at the fourth and

365

fifth supersaturation level in IPA the two curves intersect with each other, or in EtOAc the

366

third, fourth, fifth and sixth curves are almost indistinguishable in certain detection time

367

ranges. Thus summarising, the experimental data themselves seem to exhibit some intrinsic

368

inaccuracy.

369

Let us now evaluate the quality of the fitting of the model description of the experimental

370

curves. First of all, the comparison between experiments and simulations highlights the

371

inconsistency of some of the measurements, e.g. those at the lowest supersaturations in IPA

372

and EtOAc. It is also striking that the error values reported in Table 2 are the worst for

373

the experimental data that appear more problematic on visual inspection (see the discussion

374

in the previous paragraph). It is worth underlining that when the calculated cumulative

375

distribution and the measured one differ, this might be due to the wrong choice of the

376

statistics (see for instance the case of the highest supersaturation in IPA). However when

377

the distance between experimental curves obtained at different supersaturation is completely

378

different from what obtained with the model, this might be due to a wrong dependence on

379

supersaturation. While the models are able to describe the supersaturation dependence of the

380

measurements in ACN, they are not able to do so in all cases for the experiments in the other

381

two solvents. In each solvent the quality of the fitting differs at different supersaturation,

382

without exhibiting any clear trend, e.g. that lower supersaturation leads to larger errors.

383

The quality of the data in ACN appears, both on visual inspection and upon analysis of the

384

errors, to be better than that of those in EtOAc, which is in turn better than that of those

385

in IPA. It is however hard to judge whether the difficulty we experience in interpreting the

386

data in IPA and EtOAc is due to the data or to the model.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 38

Table 2: Nonlinear least-squares method. Operating conditions, total error for each individual series (from Eq.(29)) and total error for each solvent (from Eq.(30)). Series k 1 2 3 4 5 6 7 8 Es

ACN S Ek 1.05 0.83 1.07 1.21 1.09 2.08 1.11 0.96 1.14 2.04 1.17 1.51 0.62

IPA S Ek 1.07 4.78 1.09 2.95 1.10 2.69 1.12 3.81 1.14 3.87 1.18 1.41 1.22 3.19 1.26 4.62 1.26

EtOAc S Ek 1.07 3.52 1.08 1.27 1.10 4.59 1.11 1.57 1.14 1.45 1.17 2.55 1.18 1.51 1.00

387

4.2

Statistical analysis of nonlinear least-squares regression

388

Before introducing the second method, it is useful to evaluate the NLSM-results in light

389

of the representativity issue outlined in Section 2.2. Indeed, the least-square estimation is

390

based on an error minimisation where experimental measurements are implicitly assumed

391

to be affected only by Gaussian error (if the error were negligible, the measurements would

392

be accurate), whereas the measured empirical cumulative distribution function intrinsically

393

deviates from the underlying distribution, as clearly stated by the DKW inequality. Figure

394

2 illustrates the issue. For each solvent two series of experiments at two different super-

395

saturation levels are shown, by plotting the empirical cumulative distribution function, F

396

(symbols), and the corresponding distribution calculated with the model, P (solid black

397

line). The two series are chosen as those corresponding to the smallest (upper part of the

398

figure) and to the largest (lower part of the figure) model error (see Table 2). Besides,

399

in each diagram of Figure 2 three uncertainty regions around the experimental points are

400

shaded, namely those obtained by applying the DKW inequality for the values of P equal to

401

0.10 (the most colourful and smallest region), 0.70 (the intermediate region), and 0.99 (the

402

largest grey region). One sees immediately that the distance between the experimental and

403

the calculated distributions differs in the six diagrams, and that the calculated distributions 20 ACS Paragon Plus Environment

Page 21 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


404

are contained in different uncertainty regions. In the three upper diagrams the calculated

405

solid line is contained in the smallest uncertainty region, i.e. that corresponding to P=0.10.

406

In the three lower diagrams the calculated distribution is (barely) contained in the smal-

407

lest uncertainty region only in the case of ACN, whereas in the case of EtOAc it reaches

408

the largest uncertainty region corresponding to P=0.99 and in the case of IPA it goes even

409

beyond it.

410

Although intuitively this suggests that the quality of the model is much better in the

411

case of ACN than in the other two, could the deviation between experiments and model in

412

the case of IPA (lower diagram) simply be due to the unavoidable stochastic effects? The

413

analysis presented in the next section looks at such question from a different perspective,

414

namely by developing a method for parameter estimation alternative to the NLSM discussed

415

in Section 4.1.

416



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2: Comparison between the model (solid black lines) and the experimental points (grey symbols), accounting for the uncertainty related to the estimation of F (green areas). We have selected the dataset 4 and 5, 6 and 5, and 4 and 3 for ACN, IPA, and EtOAc, respectively, since they represent, for each solvent, the best (upper figure) and the worst (lower figures) fittings. The uncertainty regions are computed for increasing values of P, namely 0.1, 0.7, and 0.99 (from red, blue, and green to grey colour, for the three solvents). The model fitting (black solid lines) lies within the boundaries given by Eq. (9), even though, for example for the dataset in IPA, part of the calculated curve still lies outside the P = 0.99 boundaries. 417

4.3

418

Based on the considerations in the previous section, we are encouraged to propose a different

419

approach towards parameter estimation. This shall be based on the DKW inequality and its

420

underlying statistical model, and will exploit the p-value concept, with which one measures

421

the probability that a certain assumption - in our case the consistency between model and

422

measurements - is not fulfilled. This new approach towards parameter estimation is based

423

on the following statistical analysis.

424

Parameter estimation through p-value maximisation

Let us consider first, for each solvent independently, the complete set of experimental


Page 22 of 38

Page 23 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


425

series at different supersaturation levels. We intend to carry out a two-step statistical ana-

426

lysis, where we check first the ability of the model to describe each individual experimental

427

series k belonging to a specific solvent, and then the ability of the model to describe the

428

whole set of experimental series belonging to that solvent.

429

For a given set of model parameters that have been estimated, for every experimental

430

series k, one can compare the empirical cumulative distribution function and the corres-

431

ponding distribution obtained with the model by calculating the Chebyschev norm of their

432

difference, Dk , using Eq.(31).

433

The first step of the statistical analysis is that of testing a first level null hypothesis, i.e.

434

at the level of the individual experimental series, which states that the empirical cumulat-

435

ive distribution function is sampled from the corresponding distribution obtained with the

436

model. To this aim we determine a p-value, pk , for each series (consisting of Mk experimental

437

points). The p-value measures in this case the probability, under the statistical model asso-

438

ciated to the DKW inequality, that the Chebyschev norm of the difference above, Dk , would

439

be equal to or larger than its observed value 40 . Such p-value is the complement to one of

440

the probability defined by Eq. (7) and is given by:

pk = 2 exp −2Mk Dk2

(32)

p ln 2/(2Mk ), we set pk = 1, to avoid a p-value larger than one. The

441

Note that when Dk

Stochasticity in Primary Nucleation: Measuring and ... - ACS Publications

Recommend Documents