5216
J . Phys. Chem. 1989, 93, 5216-5224
Stochastic Theory of Nonequilibrium Phase Transitions. Applicatlon to Phospholipid Bilayer Membranes Istvan P. Sugar Departments of Biomathematical Sciences and Physiology & Biophysics, The Mount Sinai Medical Center, New York. New York 10029, and Institute of Biophysics, Semmelweis Medical University, Budapest, Hungary (Received: August 1I , 1988; In Final Form: February 13, 1989)
I n this paper the kinetics of phase transitions of the first kind are studied in separable systems. Separable systems show short-range order but do not possess long-range order such as one-component, one- or two-dimensional solids or multicomponent solids of any dimension. The equilibrium and nonequilibrium phase transitions of the whole separable system have been characterized by the processes that take place within an independent subsystem of the separable system. The phase transition has been considered as a stochastic process where the order parameter of the transition is the stochastic variable. The nonlinear, coupled stochastic equations have been solved numerically by means of Adam’s method. The processes were generated by three different types of temperature changes: linear heating/wling, temperature jumps, and sinusoidal temperature oscillations. It was possible to carry out a comparative study of these processes within the framework of the same model. Processes from unstable to stable states start with increasing velocity, but after reaching an inflection point at time til, the processes slow down exponentially. When the initial unstable state is close to the equilibrium, 2,’ is negligibly small compared to the time scale of the exponential decay and the whole process can be approximated by an exponential decay. This is the regime of the linear Onsager theory. Processes from metastable to stable states exhibit more complicated kinetics. The velocity of the process increases until the first inflection point at time t i l , and then the process slows down to a locally minimum velocity at t I Z .Then the velocity increases again till the third inflection point at tI3and finally slows down exponentially to zero velocity. Recently, these theoretically predicted kinetics were obtained experimentally, too, for the crystalline to gel phase transition of dipalmitoylphosphatidylcholine (DPPC) membranes. By means of the model, additional aspects of the phase transition kinetics have been studied, too, such as hysteresis, nucleation, homogeneous growth, relaxation processes, and phase and amplitude of steady-state structural oscillations. The model results are compared with the available experimental data taken from DPPC and dimiristoylphosphatidylcholine(DMPC) bilayer membranes.
1. Introduction The theory of nonequilibrium phase transitions is an extremely active area in theoretical physics. Its results are applied in a wide range of natural sciences and also in social sciences.’-’ The nonequilibrium phase transition is usually considered as a stochastic process. The respective stochastic equations are solved either by numerical integration or by different approximations applied in order to simplify the mathematical problem, for example, adiabatic elimination techniques’ and the method of partial distribution functions.6 The processes are directly simulated by the Monte Carlo method.’ There are special cases when the general or particular solutions (e.g., first passage time and its variation) of the stochastic equations are known.* Recently, Kanehisa and Tsong9 and Adamlo have developed models of the phase transition kinetics of phospholipid membranes. Kanehisa and Tsong adapted Fisher’s model” of drop condensation in order to model the whole phase transition process. However, NagleIz emphasized that one can describe only the beginning of ( 1 ) Haken, H. Synergetics. An Introduction; Springer-Verlag: New York, 1983. (2) Binder, K.; Stauffer, D. Statistical theory of nucleation, condensation and coagulation. Adu. Phys. 1976, 25, 343-396. (3) Kitahara, K.; Metiu, H.; Ross, J . Stochastic theory of cluster growth in homogeneous nuclei. J . Chem. Phys. 1975, 63, 3156-3160. (4) Langer, J. S.;Schwartz, A. J. Kinetics of nucleation in near-critical fluids. J . Rev. A 1980, 21, 948-958. (5) Gunton. J. D. The dynamics of random interfaces in phase transitions. J . Stat. Phys. 1984, 34, 1019-1037. (6) Volkenstein, M. V.; Gotlib, Y. Y.; Ptitsyn, 0. B. The kinetics of cooperative processes. Sou. Phys. Solid Srate 1962, 3, 306. (7) Binder, K. Monte Carlo Methods in Statistical Physics; SpringerVerlag: Berlin-Heidelberg-New York, 1979. (8) Goel, N. S.;Richter-Dyn, N. Stochastic Models in Biology; Academic Press: New York, 1974. (9) Kanehisa, M. 1.; Tsong, T. Y. Cluster model of lipid phase transitions with application to passive permeation of molecules and structure relaxation in lipid bilayers. J . A m . Chem. SOC.1978, 100, 424-432. ( I O ) Adam, G .Cooperative transitions in biological membranes. Proceedings of the Symposium on Synergetics; Haken, H., Ed.; Teubner Press: Stuttgart, 1973; pp 220-231. ( I 1) Fisher, M. E. Theory of condensation and the critical point. Physics 1967, 3, 255-283. (12) Nagle, J. F Theory of lipid bilayer phase transitions. Annu. Rel;. Phys. Chem. 1980, 31, 157-195.
0022-3654/89/2093-52l6$01.50/0
the nucleation process of the gel clusters by Fisher’s drop model. The Adam modello is a kinetic king model solved by the method of partial distribution functions. This model faces a paradox. If, in accordance with the experiments,l2 the model represents a phase transition of the first kind, it results in complex values for the relaxation times within an interval around the half-completion of the phase transition. On the other hand, one can fit the model to the measured relaxation time, as Adam did, but in this case the model does not give a transition of the first kind. In this paper the above problems are resolved by considering the membranes as separable systems. These comprise real systems such as pure one- and two-dimensional solids or most of the multicomponent solids which show short-range order but do not possess long-range order. The equilibrium phase transition of separable systems has been discussed in our previous paperI4 with applications to one- and two-component phospholipid membranes. Its basic results on cooperative unit size are summarized and further developed for the nonequilibrium case in section 2 of this paper. The stochastic equations of the transition process are derived in section 3. Finally, in section 4 the numerical solutions of these equations are presented for three different temperature excitations of the system: (i) linear heating/cooling, (ii) temperature jumps, and (iii) sinusoidal temperature changes. By means of these solutions one may study different aspects of the phase transition kinetics, e.g., hysteresis, nucleation, homogeneous growth, relaxation processes, and phase and amplitude of steady-state structural oscillations. The model results are compared with the available experimental data taken from the aqueous dispersion of phospholipid bilayer membranes. 2. Equilibrium Phase Transition in Separable Systems In this section the equilibrium phase transition properties of separable and nonseparable systems will be summarized. According to Alexandrowicz,13the macroscopic properties of a real (1 3) Alexandrowicz, Z. Analytical treatment of excluded volume, 11. Quasi random-flight behavior. J . Chem. Phys. 1967, 46, 3800-3810. (14) Sugar, I. P. Cwperativity and classification of phase transitions. Application to one- and two-component phospholipid membranes. J . Phys. Chem. 1987, 91, 95-1 01.
0 1989 American Chemical Society
The Journal of Physical Chemistry, Vol. 93, No. 13, 1989 5217
Stochastic Theory of Nonequilibrium Phase Transitions system of ordered or semiordered interacting units are equivalent to a substitutive system consisting of independent subsystems. The semimicroscopicbehavior of these subsystems should depend solely on the average state of the real system and should be independent of the state of the neighbor subsystems. This substitution is applicable for any real system that does not possess long-range order.13.36*37 Further, these real systems are called separable systems. In general, thermal motion of longwavelength phonons destroys the long-range order in the one- or two-dimensional solid, and the Bragg peaks of the diffraction pattern formed by the system are broad instead of sharp.38" The lack of long-range order during the phase transition in two-component solids has been shown by Kaufman and O n ~ a g e r . ~ ~ The model presented in this paper refers to the kinetics of phase transitions of the first kind in separable systems. At the phase transitions of the first kindIs there are two local minima of the specific Gibbs energy function g(S,p,T), where S is the order parameter of the transition, p is the pressure, and T i s the temperature. This is a common characteristic of the Gibbs free energy function whether the transition takes place in a separable or nonseparable system. However, in the case of nonseparable systems there is an additional characteristic: the ensamble average of the Gibbs energy is not analytic at the transition. This class of the phase transitions of the first kind is called first-order phase transition^.^^ Since in a separable system the Gibbs energy is analytic at the transition,14the transition is first kind but not first order. Concerning the general properties of the separable systems, it was pointed that at a given average state of the separable system, there is a lower bound to the independent subsystem size. Below this minimum size the subsystems cannot be considered to be independent from each other. _Consequently, the partition function of the substitutive system, Q, can be factorized into the partition functions of the smallest independent subsystems, Qsub
where fi and N are the numbers of interacting units (molecules) in the total system and in the smallest independent subsystem, respectively. The partition function of the smallest independent subsystem is
I
'I
I
I
314 3 314 5 3147
.I
, I
I
1
I
0
T[Kl
I
,
,
05
e
1
I
I)
I0
Figure 1. Size of the smallest independent subsystem for the gel-liquid crystalline phase transition of DPPC dispersions. Subsystem size N as a function of the (a) temperature T(K) and (b) fractional completion 8. These curves were calculated by means of eq 6 , A4, and A5 and by using the following experimental data? T d 5 = 314.49 K, T,,, = 314.55 K, T, = 314.525 K (see definitions in Appendix A).
Substituting eq 3 and 4 into eq 5 gives the size of the smallest independent subsystem as N = In [(I - s ) / e l / [ ( ( g ) l - (g),)/kTl
(6)
The temperature T, at the half-completion of the phase is called the midpoint temperature. Note, transition, Le., at 0 = the temperature Th at the maximum of the excess heat capacity curve is different from T, when the excess heat capacity curve is asymmetric (Appendix A). In theoretical papers T, is called the phase transition temperature, while in experimental papers either T,,, or Th is called the phase transition temperature. In nonseparable systems the phase transition temperature of a first-order phase transition is defined by the position of the singularity of the average Gibbs energy. At this temperature the Gibbs energies of the coexisting phases are the same. This is the Ehrenfest definition of the phase transition temperature. There is a strong correspondence between the midpoint temperature and the Ehrenfest definition of the phase-transition temperature beimplies (g), = (g)](see eq 3-5). cause e = In the middle of the phase transition, both the numerator and denominator in eq 6 become zero. The temperature derivatives of the numerator and denominator (L'Hospital rule) give N at the midpoint temperature of the transition T,: 4kTm2(de/dT)Tm
where g(S,p,T) is the specific Gibbs energy function of the real system, S is the order parameter characterizing the phase transition, p is the pressure, and T i s the temperature. We note that the effects of impurities and boundaries of the real system are included in g(S,p,T). In the phase transition range, there are two local minima of the specific Gibbs energy function at S, and S, (see, e.g., Figure 2). States around these minima are the g and I states of the system, respectively. Qfuband Q]sub contain integrals over the g and 1 states, respectively. The g and 1 states are separated from each other by state S, belonging to the local maximum of g(S,p,T) between S, and SI. Note that S,, SI,and S, depend on the actual pressure and temperature. Now the average specific Gibbs energy of g and 1 states ((g), and (g)J can be defined as
The proportion of the system in 1 phase is the fractional completion of the phase transition 8. Since the subsystems are independent from each other, the fractional completion of the phase transition for a subsystem is equal to the fractional completion of the phase transition for the whole system, 0. Thus
~
( I 5 ) Landau, L. D.; Lifshitz, E. M Course of Theoretical Physics; Pergamon: New York, 1980; Vol. 5 , Chapter XIV.
N,,, =
( ( h ) ,- ( h ) , )
(7)
where ( h ) , and ( h ) , are the average specific enthalpies of the 1 and g states, respectively. Considering that eq 7 is equivalent with the van't Hoff equation, which defines the cooperative unit size of a transition, consequently N , is equal with the cooperative unit size.14 At temperatures close to T, the average specific Gibbs energy difference can be obtained by the equation (derived in ref 14) (g)l- (g), - A L I ( T - Tm)/Tm (8) where Ahal is the molar calorimetric enthalpy change during the phase transition. Using eq 6-8 and the T,, Ahal, and e( T ) data from equilibrium calorimetry, we can determine the size of the smallest independent subsystem either as a function of temperature or as a function of the fractional completion. By assuming that the excess heat capacity curve is an asymmetric Lorentzian, analytic forms of N( T ) and N ( 0 ) are derived in Appendix A. In Figure l a and Figure 1b N is plotted as a function of temperature and fractional completion, respectively, for the case of the gel to liquid crystalline transition of DPPC membranes. Here we used the calorimetric data of Albon and Sturtevant.21 Note that because of the approximation in eq 8, N(T) is valid only at a temperature interval where I ( T - T,)/TJ 0. The spinodal temperatures have been defined in section 3 . 1 . The order parameter S , belonging to the maximum between the minima separates the order parameters into two intervals. The states belonging to the gel phase are in the 0 < S < S, interval, and the respective minimum of the Gibbs energy function is the gel minimum. The states belonging to the liquid crystalline phase are situated in the S, < S < 1 interval, and the respective minimum is the liquid crystalline minimum. The states at the global minimum represent the stable phase of the system while states around the local minimum characterize the metastable phase. At the spinodal temperature where the local minimum disappears, the metastable phase becomes unstable. By substitution of eq 6 and 16 into eq 14 and 15, the transition probabilities of the Master equations of the DPPC phase transition kinetics (eq 1 3 ) can be obtained. In the next section the numerical solutions of these Master equations are presented for different types of temperature changes, while in Appendix C the details of the applied numerical methods can be found. 4. Results and Discussion Three types of temperature-induced phase transition processes have been considered:
T ( t ) = To u
+ ut
(17)
> 0 represents linear heating, u < 0 represents cooling, u is the
scanning rate, and To is the starting temperature: T(t 0 is the temperature jump upward, and AT perature jump downward: T ( t ) = To
+ TAsin [ 2 d / t p ]
(18)
< 0 the tem(19)
sinusoidal temperature change; t p is the period time, and T, is the temperature amplitude. These temperature changes are characteristic of differential scanning calorimetric (DSC) and real-time X-ray diffraction measurement^'^-^^ of T-jump and quenching e x p e r i m e n t ~ ~ ~ - ~ ’ ~ ~ ~ (19) Caffrey, M.; Bilderback, D. H. Kinetics of the main phase transition of hydrated lecithin monitored by real-time x-ray diffraction. Biophys. J. 1984, 45, 627-63 1. (20) Cunningham, B. A,; T.-Lis, W.; Quinn, P. J. Time resolved x-ray diffraction measurements of phosphatidylcholine phosphatidylethanolamine mixtures: effect of phosphatidylethanolamine acyl chain length. J . Colloid Interface Sci. 1988, 1 2 1 , 193-197.
5220
The Journal of Physical Chemistry, Vol. 93, No. 13, 1989
Sugar
and of multifrequency calorimetric measurements.28 We should also mention the volume perturbation calorimetric measurements where the system is excited by sinusoidal volume ~ h a n g e . ~ ~ , ~ ~ 4 . 1 . Linear HeatinglCooling. In Figure 3a the calculated fractional completions of the phase transition 0 (defined by eq C I ) are plotted as a function of the actual temperature at different scanning rates u . Within the range of the considered scanning rates, the hysteresis is proportional to the logarithm of scanning rate (Figure 3b). The hysteresis Hys is defined by Hys(u) = T(0 =
X,U) - T ( & ] / ~ , - u )
(20)
where T(O=I/,,u) is the temperature at the half-completion of the phase transition when the heating rate is u. This nonequilibrium midpoint temperature depends on the scanning rate. T ( C ~ = ~ / ~is, the - U nonequilibrium ) midpoint temperature in the cooling mode. Since at a zero scanning rate the hysteresis is finite, the found logarithmic relationship should fail at small scanning rates. Recently, Yang and Nagle33measured the crystalline to gel phase transition of DPPC membranes at different scanning rates. The shift of the midpoint temperature from the equilibrium T,,, (=268.8 K) is proportional to the logarithm of the scanning rate except for the slowest scanning rate (Figure 3b). The fitted equations to the linear section of the measured and calculated hysteresis curves are
In
u
= 0.28IHys - 275.951
(21 1
for the measured dataj3 of the crystalline to gel phase transition
In
u
= 5.9[Hys - 276.311
a 'Ol-
(22)
for the calculated gel to liquid crystalline transition, where the hysteresis is measured in kelvin and the prefactors are given in K-I. Equations 21 and 22 imply that u(Hys l)/u(Hys) = 365 for the gel to liquid crystalline transition and u(Hys + l)/a(Hys) = 1.3 for the crystalline to gel phase transition. Thus 1 K increase in the hysteresis requires a 276 times faster scanning rate for the gel-liquid crystalline phase transition than for the crystalline-gel transition. We propose now a general relationship between the scanning rate and the hysteresis which is also valid for the case of slow scanning rate. In the case of separable systems the hysteresis
+
(21) Albon, N.; Sturtevant, J. M. The nature of the gel to liquid crystalline transition of synthetic phosphatidylcholines. Proc. Natl. Acad. Sci. U . S . A . 1978, 75,2258-2260. (22) Eigen, M.; De Maeyer, L. C. Tech. Org. Chem. 1963,8,845-1054. (23) Trauble, H. Phasenumwandlungen in Lipiden, Mogliche Schaltprozesse in biologischen Membranen. Naturwissenschafren 1971, 58, 277. (24) Tsong, T . Y . ; Kanehisa, M. I. Relaxation phenomena in aqueous dispersion of synthetic lecithins. Biochemistry 1977, 16, 2674-2680. (25) Gruenewald, B.; Blume, A.; Watanabe, F. Kinetic investigations of the phase transition of phospholipid bilayers. Biochim. Biophys. Acta 1980, 597, 41-52. (26), Genz, A,; Holzwarth, J. F. Dynamic fluorescence measurements on the main phase transition of dipalmitoylphosphatidylcholine vescles. Eur. Biophys. J 1986, 13, 323-330. (27) Teissie, J. Fluorescence temperature jump relaxations of dansylphosphatidylethanolamine in aqueous dispersions of dipalmitoylphosphatidylcholine during the gel to liquid crystalline transition. Biochim. Biophys. Acta 1979, 555, 553-551. (28) Mayorga, 0. L.; van Osdol, W. W.; Lacomba, J. L.; Freire, E. Frequency spectrum of enthalpy fluctuations associated with macromolecular transitions. Proc. Natl. Acad. Sci. U.S.A. 1988, 85, 9514-9518. (29) Johnson, M. L.; van Osdol, W. W.; Biltonen, R. L. The measurement of the kinetics of lipid phase transitions: a volume perturbation kinetic calorimeter. Merhods Enzymol. 1986, 130, 534-551. (30) van Osdol, W. W. Kinetics of the main phase transition in phospholipid vesicles, Ph.D. Thesis; University of Virginia, Charlottesville, 1988. (31) Nagle, J. F.; Scott, H. L. Lateral compressibility of lipid mono- and bilayers. Theory of membrane permeability. Biochim. Biophys. Acta 1978, 513, 236-243. (32) Doniach, S. Thermodynamic fluctuations in phospholipid bilayers. J . Chem. Phys. 1978, 68, 4912-4916. (33) Yang, C. P.; Nagle, J. F. Phase transformation in lipid follow classical kinetics with small fractional dimensionalities. Phys. Reu. A: Gen. Phys. 1988, 37. 3993-4000.
315
314
316
T[Kl
Figure 4. Linear heating and cooling processes. In these calculations the same model parameters were used as in Figure 3a. (a) The average specific Gibbs energy ( g ) as a function of the actual temperature T. (b) The fluctuation of the disorder parameter ( A S ) a s a function of the actual temperature. Solid line: the scanning rate u = +10"(K), dashed line: u = -10"(K), where u = u,,,/Y.
vanishes at zero scanning rate, Le., Hyso = 0 (see Figure 3a). During this infinitely slow transition the order parameter possesses equilibrium (Boltzmann) distribution, and consequently the transition is an equilibrium transition. In nonseparable systems the hysteresis does not vanish. Even in the case of a very slow scanning rate the hysteresis approaches the difference of the spinodal temperatures, Le., Hyso = T,, - T,, > 0. This is a nonequilibrium phase transition at zero scanning rate. The following formula results in the proper hysteresis at slow scanning rate and satisfies the previously found logarithmic rule at faster scanning rates: u = i(exp[B(Hys)] - exp[B(Hyso)])
(23)
where and B are model parameters. In Figure 4a the ensemble average of the specific Gibbs energy ( g ) is plotted as a function of the actual temperature at heating and cooling modes. (g) is defined by
where T ( t ) is given by eq 17. In thermodynamics the intercept of the specific Gibbs energy curves obtained by quasi-stationary heating and quasi-stationary cooling determines the transition temperature of first-order phase transitions. In separable systems this temperature corresponds to the equilibrium midpoint temperature of the phase transition of the first kind (see section 2). In our case the intercept of the calculated Gibbs energy curves is at 314.95 K (Figure 4a), and really, this temperature is equal with the midpoint temperature of the phase transition at zero scanning rate (see Figure 3a). Above this phase transition temperature in Figure 4a the solid line represents the specific Gibbs energy of the metastable gel phase, while at lower temperatures the dashed line gives the specific Gibbs energy of the metastable liquid crystalline phase. In the metastable phases, the larger the deviation from the phase transition temperature, the faster the transition into the stable phase. Figure 4b shows the fluctuation of the order parameter ( A S ) as a function of the actual temperature, where M n=O
and M
Comparing Figure 4a and Figure 4b, the fluctuations are enhanced close to the phase transition temperature, especially in the metastable phases. In the case of lipid membranes the functional importance of the large fluctuations around the phase transition
Stochastic Theory of Nonequilibrium Phase Transitions log
Ii dB -
I
2,,2
The Journal of Physical Chemistry, Vol. 93, No. 13, 1989 5221
N
900 700
dT
500 300
I
.
.
,
8
,
.
,
,
314
#
.
,
.
.
_
314 5 T
[Yl
Figure 6. Characteristic times of the quenching processes. Starting temperature: To = 3 16 K; T(K) is the final temperature of the quenching process. Line A: the place of the first inflection point til of O ( t ) kinetics as a function of the final temperature T. Line B, C: the place of the second and third inflection points of O ( t ) kinetics, ti2and ti3, respectively. Dashed line: half-time of the final exponential decay of the quenching process as a function of the final temperature T. Figure 5. Combined plot of relaxation time, subsystem size, and d6/dT curves. Solid line: the equilibrium dO/dT curve. Dashed line: equilibrium curve of subsystem size N(7‘). Dash-dotted line: half-time of the relaxation processes 71/2 as a function of the starting temperature To, the half-time is scaled by the characteristic frequency (7112 = U 7 1 / 2 ( r a l ) ) . The temperature jump is AT = +0.01 K. (a) The model parameters were taken as in Figure 1. (b) Hypothetical model parameters: six times broader transition than in the case of (a): T,, = 314.7625 K, T+o,5= 3 15.1225 K, Th = 3 14.9725 K. The asymmetry A was kept constant (see Appendix A). temperature was discussed in connection with the enhancement of ion permeability of the membranes at Tm31332 4.2. Small Temperature Jumps. In these calculations the temperature jump was kept constant ( A T = fO.O1 K) and the starting temperature To has been varied. In accordance to the linear Onsager theory of processes, these small perturbations decay exponentially (not shown). In Figure Sa the half-time of the relaxation process T l / 2 is plotted against the starting temperature To. The half-time of an exponential decay T1/2 is proportional to the relaxation time of the process 7, where T I / 2 = T In 2. In agreement with the experimental data23-27the relaxation time has a maximum at Th. rh is the temperature at the maximum of the excess heat capacity curve. The relaxation processes slow down a t Th because the fluctuations are maximal at this temperature. However, the calculated maximum of the relaxation time curve is much sharper than the experimental ones. In Figure 5a, T l / 2 spans 5 orders of magnitudes while in the dynamic fluorescence measurements of Genz et a1.26it changes by only 2 orders of magnitudes. But this is understandable because our calculation was carried out on the basis of Albon and Sturtevant’s data2 when we constructed the N(B) curve (see Appendix A). Since Albon and Sturtevant’s sample was at least 99.94% pure, the temperature range of the transition was very narrow (see solid line in Figure Sa). The relaxation time measurements have been performed on less pure samples containing fluorescent probe molecules, and therefore the temperature range of the transition is broader. In ordFr to illustrate the effect of impurities on the relaxation time curves, in Figure 5b the half-width of the heat capacity curve was set six times larger than in the case of Albon and Sturtevant’s measurements, and now the relaxation time values calculated by means of this hypothetical heat capacity curve span only 1 order of magnitude. 4.3. Large Temperature Jumps. In this set of calculations the starting temperature To has been set above the equilibrium midpoint temperature T , of the phase transition ( T o = 316 K,
314 9
04
314 89
2 109
I 09 I 0L 314 96
i!b
8 04 02
500
I000
t
1500
2000
Figure 7. Processes generated by large temperature jump. Starting temperature To = 316 K; final temperature T(K) is shown at each curve. Real time is scaled by the characteristic frequency Y. The model parameters were taken as in Figure l . (a, b) Relaxation processes and nucleation processes. (c) Homogeneous growth processes. (d) The final exponential part of the nucleation process; enlarged part of the O(t) curve at T = 314.7 K.
T , = 314.95 K) and AT has been changed systematically. Further, only the quenching processes will be considered since the upward temperature jump processes show qualitatively the same features. The beginning of the quenching processes does not depend too much on the final temperature. In each case the velocity of the process increases up to a certain time til. In Figure 6, line A shows t,] as a function of the final temperature T. From time til the quenching process starts to slow down. Thus from a mathematical point of view, the process has an inflection point at t , , . Depending on the final temperature T , three qualitatively different processes take place: (i) relaxation processes, (ii) nu-
5222
The Journal of Physical Chemistry, Vol. 93, No. 13, I989
cleation processes, and (iii) processes of homogeneous growth. Figure 7 shows the time courses of these quenching processes at different final temperatures T . 4.3.1. Relaxation Processes. When the final temperature T > T , ( ~ 3 1 4 . 9 5K ) , the quenching process is very close to exponential decay, as in the case of small temperature jumps. The initial accelerating part of the relaxation process is negligibly short relative to the time'scale of the whole process. The relaxation processes do not lead to phase transition since the final temperature is still higher than the phase transition temperature, but during the processes the proportion of the gel phase increases. When the final temperature Tis close to T , (=314.95 K), the relaxation process starts to deviate seriously from the single-exponential decay. At 3 14.95 K a long linear part appears in the time course of the process. In Figure 6 at T = 3 14.95 K, a vertical line represents the time interval where the time course is linear. Since dZB/dt2is zero at this linear part of the curve, one can consider these points of the kinetics as a continuous series of inflection points. Otherwise, this is the limit process between the relaxation and nucleation processes. 4.3.2. Nucleation Processes. When the final temperature is between T , and 3 14.3 K, the B(t) kinetics have a quite complicated character. Up to til the velocity of the process increases. Then the process slows down to a locally minimal velocity at time ti*. At this point 0 > 0.5; i.e., the majority of the system is still in the liquid crystalline phase (see Figure 7b). From ti2 the velocity increases again until the third inflection point of the process is at ti3, and finally the process slows down exponentially to zero velocity (see Figure 7d). The phase transition process has been completed. In Figure 6 lines B and C show tiz and ti3, respectively, as a function of the final temperature T. The dashed line in Figure 6 is the half-time of the final exponential decay of the quenching processes (again as a function of T). When the half-time of the final exponential decay is much smaller than ti3, the appearance of the final equilibrium state looks instantaneous. This is the typical feature of the kinetics in Figure 7b. Figure 7d is an enlargement of the final relatively fast exponential decay at 314.7 K. But at T = T,, the half-time of the final exponential decay is comparable with ti3 and we can see the final exponential decay (see Figure 7a at T = 314.89 K). These processes are called nucleation processes because transitions from metastable to stable phases take place. A latency period precedes the phase transition. This period of time (=ti3) is the lifetime of the metastable phase, while the states of the subsystems are distributed around the local minimum of the Gibbs energy function rather than around the global minimum. When, as a result of the fluctuations, the critical value of the order parameter S , (defined in section 2 ) is reached, the transition may take place in the considered subsystem. At the lower limit temperature of the nucleation processes ( T = 3 14.3 K ) , again a long linear part appears in the time course of the process as in the case of the upper limit temperature of the nucleation processes ( T = 314.9 K). This is the limit process between the nucleation processes and homogeneous growth processes. Note that the lower limit temperature of the nucleation process is higher than the lower spinodal temperature T,, (=314 K). Thus, the quenching process loses the characteristics of the nucleation processes albeit the local minimum of the Gibbs energy function belonging to the liquid crystalline phase still exists. This is the case because the local minimum is very shallow. Its depth is comparable with the thermal energy unit kT and it cannot keep the system in the metastable state. 4.3.3. Homogeneous Growth. When the temperature drops below 314.3 K. the O ( t ) curves have only one inflection point: at til (Figure 6). Since til is comparable with the time scale of the whole process, we can see the initial accelerating part of the quenching process (Figure 7c). The liquid crystalline to gel transition is fast relative to the nucleation-type transition process. The liquid crystalline phase is unstable from the very beginning, a n d the homogeneous growth process is driven primarily by
Sugar (Si
Figure 8. Processes generated by large temperature jump. The parameters of the calculationscorrespond to that of Figures 1 and 7. The upper curves and the lower curves represent the time dependence of the average order parameter (S') and that of the fluctuation of the order ( A S ) , respectively. The final temperatures T ( K ) are shown at each curve. (a, b) Relaxation processes and nucleation processes. (c) Homogeneous growth processes.
deterministic factors Fd rather than by random factors F, (see eq 9). 4.3.4. Experimental Evidence. Can we measure the above kinetics of quenching processes, especially the complicated time courses of the nucleation processes? The main experimental obstacle is the very long lifetime of the metastable phase relative to the fast phase transition process. In order to avoid this problem, one has to quench the system almost to the lower limit temperature of the nucleation processes. Here the characteristic features of the nucleation processes still exist and the lifetime of the metastable state is minimal (Figure 6). Recently, Tristram-Nagle et measured the kinetics of the gel to crystalline phase transition of DPPC membranes. The fractional completion of the transition as a function of time possessed three inflection points when the system was quenched to 278.05 K. This is below the equilibrium phase transition temperature of the gel to crystalline temperature T , (=286.95 K ) . Further lowering the final temperature, the kinetics possessed only one inflection point. The measured O(t) kinetics show the initial accelerating region of the quenching processes. These quenching experiments are in close agreement with the predictions of our model. On the other hand, the classical Avrami theory of phase transition kinetics34is not able to predict the presence of inflection points in the kinetics. For completeness Figure 8 shows the change of the average order parameter (S) and that of the fluctuation (AS) (defined by eq 25) during the quenching processes. The common feature is that phase transition processes (nucleation or homogeneous growth) are coupled with the enhancement of the fluctuations, with a maximum at about the half-completion of the transitions, while in the case of relaxation processes, when phase transitions are not involved, the fluctuations change monotonically. 4.5. Sinusoidal Temperature Oscillation. In the case of sinusoidal temperature oscillation, (see eq 19) after 10-20 periods, the response of the system will be periodic, too. For example, the (34) Avrami, M. J . Chem. Phys. 1939, 7 , 1103; Ibid. 1940, 8, 212; Ibid. 1941, 9, 177. (35) Freire, E. Personal communication. (36) De Gennes, P. G.Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca-London, 1979. (37) Fisher, M. E. Personal communication. (38) Peierls, R. E. Helu. Phys. Acta 1923, 7 , 81. (39) Landau, L. D. Phys. Z . Sowjetunion 1937, 11, 26. (40) Mermin, N. D. Crystalline order in two dimension. Phys. Rev. 1968, 176, 250-254. (41) Kaufman, 9.; Onsager, L. Phys. Reu. 1949, 76, 1244. (42) Gardiner, C. W. Handbook of Stochastic Methods; Haken, H., Ed.; Springer-Verlag: Berlin-Heidelberg-New York, 1985. (43) Tristram-Nagle, S.;Wiener, M. C.; Yang, C. P.; Nagle, J . F. Kinetics of the subtransition in dipalmitoylphosphatidylcholine. Biochemistry 1987, . . 26, 4288-4294. (44) Callen, H. 9. Thermodynamics; Wiley: New York, 1961.
The Journal of Physical Chemistry, Vol. 93, No. 13, 1989 5223
Stochastic Theory of Nonequilibrium Phase Transitions
-3
-4
-2
-1
log [ I / t p l
I:O~
0
b.
-4
L
”
-2
T[Kl
log I / t p )
Figure 9. Steady-state responses to the sinusoidal temperature oscillations. The parameters of the calculations correspond to that of Figures 1 and 7 . The amplitude of the temperature oscillation TAwas set to 0.01 K. The period time tp of the oscillation is scaled by the characteristic frequency u. (a) The logarithm of the response amplitude 8, and (b) the phase delay A 9 are plotted as a function of the logarithm of 1 / z p . The starting temperature Tois shown at each curve. (c) The response amplitude is plotted as a function of the starting temperature at different period time tp shown at each curve.
fractional completion O(t) oscillates around the equilibrium value belonging to the starting temperature 6( To):
6(f)
E
TO)
+ 6,
Sin
[ [ 2 d / t p ] -‘Ap]
(26)
where the amplitude of the response BA and the phase delay A p of the oscillation are complicated functions of the starting temperature To,the amplitude of the temperature oscillation TA,and the period of the temperature oscillation tp. In the Figure 9 these functions are shown where TA is 0.01 K. In Figure 9a the amplitude of the response 6, is plotted as a function of the excitation frequency 1 / t p at different starting temperatures To. The amplitude decreases with increasing frequency and with To approaching Th as well. (Th is defined by the maximum of the equilibrium excess heat capacity curve in Appendix A.) The response amplitude becomes insensitive to the excitation frequency above a critical frequency, uc (Y, is defined by the minimum of the response amplitude). In Figure 9b the phase delay Ap is plotted against the excitation frequency at different starting temperatures. The general property of these curves is the cutoff of the phase delay. The cutoff frequency and the above defined critical frequency uc coincide (cf. Figure 9a with Figure 9b). The decrease of the phase delay is very sharp at each temperature, and the cutoff frequencies are grouped within a narrow interval. The phase delay decreases also toward the lower frequencies, but the decrease is not sharp at all, and the slope depends strongly on To. In Figure 9c the response amplitude is plotted as a function of the temperature To at different excitation frequencies. As was mentioned before, the closer T is to Th, the smaller is the response amplitude. Temperature and frequency dependence of the amplitude response associated with dimiristoylphosphatidylcholine (DMPC) vesicles have been measured by a recently developed multifrequency calorimeter.28 In qualitative agreement with the prediction of our model, the response amplitudes decrease with increasing frequency and also with temperatures approaching Th,but there is no evidence yet for the predicted critical frequency uc. 5. Conclusions The presented stochastic theory has been developed for the description of the kinetics of phase transitions of the first kind
in separable systems. The transition process in the system as a whole has been characterized by the processes that take place in an independent subsystem. The model made it possible to compare relaxation processes, nucleation processes, and homogeneous growth processes together. The general methods of the theory have been applied for the transition processes of one-component phospholipid membranes. The processes have been induced by three types of temperature changes, but the model can be used for the quantitative description of nonequilibrium phase transitions induced by the change of any other thermodynamic variable. Characteristic, experimentally detected features of the phase transition kinetics of DPPC and DMPC membranes have been obtained by means of the model, such as the exponential relationship between the scanning rate and the hysteresis of the phase transition, the maximum of the structural fluctuations at the half-completion of the transition, the maximum of the relaxation time of the small perturbations at temperature Th defined by the maximum of the excess heat capacity curve, the complicated phase transition kinetics with one or three inflection points, and the minimal response amplitude at temperature Th in the case of periodic excitation.
Acknowledgment. I thank R. L. Biltonen, T. E. Thompson, M. L. Johnson, W. W. van Osdol, M. Lacker, J. Bond, I. Salgo, and E. Freire for their helpful discussions. This research was supported by the U.S. Army Research Office and Dr. T. E. Thompson’s USPHS Grant GM-I 4628.
Appendix A Derivation of an Analytic Form for the Subsystem Size. By means of eq 6 and the experimentally obtained equilibrium excess heat capacity curve Cp(T),one can determine the size of the smallest independent subsystem N as a function of the temperature T or fractional completion 8. In the case of one-component systems, Cp(T ) dB/dT is a fairly good a p p r ~ x i m a t i o n land ~ dO/dT can be approximated by an asymmetric Lorentzian:
-
X..at T
C
Th
X + at T
> Th
where X* = ( T - Th)/(T*o,s - Th); 6’h = dO/dT at T = Th; Th is the temperature at the maximum of the C ,( 7‘) curve; T4,5and T+o,sare the lower and upper temperatures at the half of the peak of the heat capacity c y v e , respectively. The asymmetry of the heat capacity curve A is characterized by
2 = (T+o.s- Th)/(Th - T-0.5)
(‘42)
According to Albon and Sturtevant’s calorimetric data2’ on DPPC membranes, T4,5= 314.49 K, T+o,s= 314.55 K, Th = 314.525 K, and A = 0.71. The slope of the fractional completion at Tt is 6’h = 10.61 K-l obtained from the following normalization condition:
By integrating eq A1 we obtain the fractional completion as a function of the temperature: O ( T ) = 6’h{[T*0,S- Thl arctan
(x*)- r[T-0.5- Thl)
(A4)
and the inverse function T(6) is
where in eq A4 and A5 the minus sign refers to T < Th and the plus sign to T > Th. Substituting eq A4 or AS into eq 6, one obtains N as a function of T or 6, respectively (see Figure 1). These maximum curves are always asymmetric to the position of the maximum even if the heat capacity curve is symmetric.
5224 The Journal of Physical Chemistry, Vol. 93, No. 13, I989 Appendix B Approximation of the Fokker-Planck Equation by a Master Equation. It is a fundamental result of the theory of stochastic processes4*that a diffusion process can always be approximated by a jump process. In our case a partial differential equation ( Fokker-Planck equation) describes the process of the phase transition, where the time t and the order parameter S are continuous variables. In order to solve a set of ordinary differential equations rather than a partial differential equation, we approximate the Fokker-Planck equation by a birth-death Master equation. In the Master equation the time t is continuous variable and the order parameter S is discrete variable. If the FokkerPlanck equation has the general form aP(S,t) -at
- - a M S ) ~ ( ~ , t+) l-I as 2
a2[B(s) p(S,t)l as2
(B1)
then the Master equation which approximates eq B1 is dPS, dt
- %*Ips*,
+ Gs+lps,l
-
[%-" + %,IPS"
(B2)
where the ws, transition probabilities are in the following relationships with A(S) and B(S):42 Gsn = [A(S,)/21
+ B(S,)/2I2]
(B3)
Gsn = [-A(S,)/21
+ B(S,)/212]
(B4)
In eq B3 and B4 1 is the distance between the discrete values of the order parameter; i.e., S, = nl. Psn(t)is the probability that the order parameter is S, at time t . Esmis the probability of an S, Sn+lprocess during dt time; Gsn is the probability of the S, S,, process during dt time. In the stochastic equation of the phase transitions (eq 12)
--
ac as
A(S) = h -,
B(S) = Q
(B5)
where G(S,p,T) is the Gibbs energy function of a subsystem; h is a proportionality factor and Q is the diffusion coefficient. By substituting eq B5 into eq B2-B4 we get the Master equation of the phase transitions. Note, if h = Q/(kT) and the step size 1 is small enough, the transition probabilities defined by eq B3 and B4 are equivalent with the following definition of the transition probability functions: GS" = v exp(-[G(S,+l,p,T) - G(S,,p,T)l/kT)
(B6)
Gs"=
037)
v
exp(-[G(S,l,p,T)
- G(S,,P,T)I/~T)
where u = Q/(212). The analogy between eq B3-B5 and eq B6 and B7 can be pointed out by using the approximation exp(-[G(S,+,,p,T) - G(S,,p,T)I) = 1 -
(B8)
These forms of the transition probabilities are widely used in Monte Carlo type simulations of the phase transitions.2 The stationary solution of the Master equation eq B2 can be obtained by mere summation if the principle of detailed balance is fulfilled. The principle of detailed balance requires that there are ab many transitions per second from state S, to state S,+, as from S,+, to S , by the inverse process, or in mathematical form &$Sn
= ~s"+l~s,+l
(B9)
Finally, the stationary solution of the Master equation eq B2 is1,42
ps, =
exP(-[G(S,A T )/ k T l )
(BIOI
Zexp(-[G(S,,p,T) / k T l )
,=O
where the summation goes through the whole set of order pa-
Sugar rameters. As required, the obtained stationary distribution is the Boltzmann distribution. Appendix C The Numerical Methods. In order to solve the highly nonlinear equations of the phase transition knetics, Adam's method has been used (IMSL Library routine: DGEAR). The precision was set to 10" in each run. In order to reduce the computation time, only the (0,0.35) interval of the order parameter S was divided into 100 equal sections and S = 0.35 had been considered as a reflection pont of the stochastic process8 rather than S = 1. We note that the most probable liquid crystalline state is at S = SI(-0.2) (SI is defined in section 2). This simplification did not alter the result since in the neglected (0.35,l) interval, the Psn(t) probabilities were less than in each considered case. Thus in the boundary conditions in section 3.2, S , = 0.35 had been used, actually. In all of the computations the time t was scaled by the characteristic frequency u and was given in dimensionless units ( t = ut,,l). The actual value of the fractional completion of the phase transition had been calculated by the following formula, which except the last term, is analogous to eq 5 : M
n=u+l
-n
where 1 = S, - S,,; S, is the order parameter at the maximum of the specific Gibbs energy function (defined in section 2); u = int (S,/l) is the number of mesh points in the (O,S,) interval. The function int truncates real to integer number, e.g., int (1.85) = 1. The role of the last term in eq C1 is the following. Without this term O(t) may possess stepwise changes. Namely, when the temperature is changing, S, changes too and the lower limit of the summation u may change by f l . The last term in eq C1 is a linear interpolation which smoothes out the stepwise changes. The determination of O ( t ) is based on the existence of the double minima of the specific Gibbs energy function. But the g minimum disappears at the upper spinodal temperature T,, and the 1 minimum disappears at the lower spinodal temperature T,, (defined in section 3.1). Can we define O ( t ) out of the region of the spinodal temperatures? At the spinodal temperature the disappearing local minimum of the Gibbs energy function and the local maximum at S , coincide, resulting in an inflection point. The order parameter belonging to this inflection point Si(p,T)should take over the role of S, in eq C 1. By means of this definition, in the case of infinitely slow heating/cooling I9( T ) is changing continuously even at the spinodal temperatures as is expected for phase transitions in separable systems (where N C m). By contrast, O ( T ) changes discontinuously at the spinodal temperatures during infinitely slow heating/cooling for nonseparable systems (where N = a), The initial conditions of the stochastic equations eq 13 are given by the equilibrium distribution of S , at the starting temperature To : PS"(0) =
exp(-[N(@ g(S,,P>To)/kTOl) M
(C2)
i=O
Since in eq C2 I9 itself depends on the Ps,(0) distribution (see eq C l ) , we followed a self-consistent procedure in order to determine the equilibrium distribution at To: (i) Choose a trial subsystem size of N , . (ii) Determine the respective fractional completion 19,by using first eq C2 and then eq C 1. (iii) Determine the second trial value of the subsystem size by means of eq 6, N(8,) = N,, and follow the procedure till self-consistency. We used the above procedure determining the fractional completion during infinitely slow heating and cooling as well. Registry No. DPPC, 2644-64-6.
