Alkalmazott potenciálmodellek és elvégzett számítások

A kalciumcsatornák pontos molekuláris szerkezete egyelőre nem ismert. Az általánosan elfogadott feltevés az, hogy kalciumcsatornák jellemzően nagy kalciumra vonatkozó szelektivitásáért egy, a transzmembrán fehérje pórusán belüli szűk térrész, az ún.

szelektív szűrő aminosav oldalláncainak négy, negatívan töltött karboxil-csoportjai a felelősek. A kalciumcsatornák szerkezetének modellezésére jelenleg az irodalomban többnyire valamilyen egyszerűsített modellt alkalmaznak. Számításaimban azt a Boda és munkatársai által kifejlesztett modellt [138-142] használtam, amelyben a szelektív szűrő lényegi részét képező a karboxil-csoportok oxigénjeit -1/2 parciális töltésű és 0,28 nm átmérőjű merevgömbök reprezentálják, amelyek a szimuláció során nem léphetnek ki a szelektív szűrűből (a szelektív szűrő z-irányú határai az oxigénekre nézve merev falakként funkcionálnak, ld. 26. ábra). A fehérje egyéb részeit merev, áthatolhatatlan falak modellezik. Az elektrolit oldatokat az előző fejezetben is említett primitív modell reprezentálja [137], ahol az oldószer implicit módon, a vízre jellemző, εr = 78,65 relatív permittivitással vehető figyelembe, a tömbfázis oldataiban és az ioncsatorna belsejében egyaránt. Ez a kontinuum-modell tehát elhanyagolja azt, hogy a csatornában mások a

polarizációs hatások, mint a tömbfázisban. A modellben az elektrolit-oldat nátrium-, kalcium- és klorid-ionjai töltéssel rendelkező merevgömbök, amelyek átmérői sorrendben

Na+

d = 0,19 nm, dCa2+= 0,198 nm és dCl−= 0,362 nm, parciális töltései pedig természetesen

Na+

qe = +1, qeCa2+= +2, illetve qeCl−= -1. A 26. ábra szemlélteti a szimulációs doboz elrendezését, melyet z-irányban merev falak határolnak. A szimulációk során a „minimum image” konvenciót z-irányban is figyelembe vettem, de a kontrollcellák és a szélső merev falak között 0,4 nm széles átmeneti tartományt definiáltam (a dinamikus cella és a kontrollcellák között nem volt átmeneti tartomány). Az rmax szimulációs paraméter meghatározása csak a csatorna szelektív szűrőjét és „előszobáját” tartalmazó dinamikus cella térrészében mért s úthosszak alapján történt, de az (56) egyenlettel való átlagképzés során a karboxil-csoportok oxigénjeit nem tekintettem külön komponensnek. A rendszer teljes Coulomb-energiáját a (8) összefüggéssel számítottam. Az elektrolit-oldatok komponenseinek kémiai potenciálját a Boda és munkatársai által kidolgozott iterációs eljáráson alapuló módszerrel számítottam [146, 147].

Az elvégezett számítások hossza 10¹⁰ szimulációs lépés volt, melynek felét a nagykanonikus kontrollcellák részecskeeltávolítási és -beillesztési kísérletei tették ki. A szimulációs doboz bal oldali kontrollcellájában a kationok összkoncentrációja mindig c_bk = 0,1 mol/dm³ volt. Az anyagtranszport hajtóerejének megfelelő elektrokémiai gradienst a jobb oldali kontrollcella összkoncentrációjának (cjk), továbbá a kontrollcellák kationarányának különböző beállításaival definiáltam. A két kontrollcellában a kationarány mindig megegyezett, de a jobb oldali kontrollcellában a kationok összkoncentrációja mindig kisebb volt a bal oldalihoz képest. A használt koncentrációarányok a következők voltak: αc = 1, αc = 1/3, αc = 1/10, αc = 0, ahol αc = cjk/cbk. (Az αc = 1 eset az egyensúlynak felel meg, míg az αc = 0 eset a legnagyobb hajtóerőt feltételezi, amikor a jobb oldali kontrollcellában végtelenül híg oldat van. Ez utóbbi megvalósítható úgy, hogy a komponensek kémiai potenciáljaihoz extrém nagy negatív értékeket rendelünk, és így a részecskebeillesztés itt gyakorlatilag soha nem lesz elfogadott.)

26. ábra. A kalciumcsatorna egyszerűsített modelljének szerkezeti felépítése.

5.2.3. Eredmények

Az ioncsatornák szelektivitását kétféleképpen lehet definiálni: (1) a csatornában mért átlagos betöltöttségek hányadosával és (2) a csatornán átfolyó anyagáramok hányadosával. A csatorna betöltöttsége (n) és a komponensek árama (Jz) természetesen nem függetlenek egymástól. Ha a csatorna egy adott térfogategységére vonatkozóan ν fejezi ki egyetlen ion egységnyi időre eső átlagos elmozdulását, azaz az egyedi ion átlagsebességét,

akkor ν_z (az egyedi ion átlagsebességének z-irányú komponense) és n (az adott ionos komponensre vonatkozó betöltöttség) szorzata adja a Jz áramot az adott ionos komponensre.

27. ábra. A csatorna betöltöttségének (n), és az ionáramok (Jz) hányadosai a móltört függvényében (x a CaCl2 móltörtje a kontrollcellákban).

A 27. ábra mutatja a szelektivitásokra kapott értékeket a CaCl₂ tömbfázisbeli (kontrollcella) móltörtjének függvényében. A legfontosabb következtetés, melyet levonhatunk az eredményekből az, hogy a csatorna affinitását jellemző, betöltöttségekkel kifejezett szelektivitás a kalcium-ionra mindig nagyobb, mint az ionáramok hányadosaként definiálható szelektivitás. Ez azt jelenti, hogy a kalcium-ion a nátrium-ionénál általában nagyobb átáramló mennyisége ellenére előszeretettel tartózkodik a csatornában, és nem szívesen hagyja el azt. Az eredmény mindenképpen érdekes, hiszen az ioncsatornák szelektivitásának mértékére eddig többnyire csak a csatorna betöltöttsége alapján tudtak következtetni. Érdemes megjegyezni, hogy míg a csatorna adott keresztmetszetére

addig a tömbfázisban mindig a nártium-ion mozgékonyabb (diffúziós állandója például az x

= 0,2 és 0,8 móltörteknél kb. 2,5-ször nagyobb, mint a kalcium-ioné). A 27. ábra eredményei alapján megállapítható továbbá az is, hogy a csatorna sugarának (R) csökkenésével és a hajtóerő növekedésével a csatorna egyre szelektívebb lesz a kalcium-ionra.

Az αc = 1 egyensúlyi esetnek megfelelő számítások betöltöttségre vonatkozó eredményeit összehasonlítottam a Boda és munkatársai által számított, hagyományos nagykanonikus MC szimulációkban kapott adatokkal és szinte tökéletes egyezést találtam.

Az iontranszport mechanizmusának mélyebb megértéséhez megvizsgáltam a csatorna belsejében számított tulajdonságoknak az áramlás irányában vett eloszlását (az eloszlások meghatározásához a csatornát z-irányban 0,05 nm széles térfogategységekre osztottam, és ezeken belül gyűjtöttem a keresett mennyiségek átlagát). A 28. ábrán láthatók a csatorna betöltöttségének (a térfogategységek számsűrűségének a z-irányra merőleges felülettel vett szorzata), az egyedi ionok z-irányú sebességének, illetve az ionáramoknak a profiljai.

Megfigyelhetjük, hogy a profilok asszimetrikusak a két kontrollcellában alkalmazott különböző összkoncentráció miatt, de ez az aszimmetria fokozatosan eltűnik, ahogy közelítünk az egyensúlyi esethez (ld. az αc = 0 és αc = 1/3 esetek közti különbség). A betöltöttségek eloszlásának elemzéséből az a fontos megállapítás vonható le, hogy a tömbfázisok összetételétől és a hajtóerő mértékétől függetlenül a szelektív szűrő közepén átlagban mindig ugyanannyi kalcium-ion tartózkodik (ld. az n(Ca²⁺) görbék maximumainak környezete). A görbe alatti területek nagyságát kiszámítva megállapítottam, hogy a csatornában átlagosan mindig 1-nél több kalcium-ion van, míg a nátrium-ionra ez az átlag legalább egy nagyságrenddel kisebb. Természetesen az n és νz mennyiségek szorzataként kapható Jz értékeinek minden egyes térfogategységben azonosaknak kell lenniük, hiszen ekkor teljesül az anyagmegmaradás törvénye (tehát a logaritmikus skálán ábrázolt n- és νz -profilok összege konstanst ad). Ez persze azt is jelenti, hogy a csatorna olyan térrészeiben, ahol az adott ionok valamiért összetömörülnek, ott z-irányban lassabban haladnak előre.

Azonban a számításokból az is kiderült, hogy ilyenkor az egyedi ionok átlagos sebessége nem csökken olyan erőteljesen, mint νz, ami azt jeleni, hogy az ionok mozgékonysága többé-kevésbé változatlan marad, ahogy az ionok a csatorna belseje felé haladnak.

28. ábra. A csatorna betöltöttségének (n), az egyedi ionok z-irányú sebességkomponensének (νz) és a komponensek áramainak (Jz) z-irányú eloszlásai (x a CaCl2 móltörtje a kontrollcellákban).

ÖSSZEFOGLALÁS

Az értekezésben bemutatott legfontosabb eredmények az alábbi tézispontokban foglalhatók össze.

Tézispontok:

A szimulációs szoftver elkészítése

1. Doktori munkám első célkitűzése az általam korábban fejlesztett C nyelven írt Monte Carlo szimulációs szoftvercsomag továbbfejlesztése volt. A platformfüggetlen, objektumorientált C++ nyelven írt új változat a kutatási céljaimnak megfelelően a különböző típusú MC lépések engedélyezésével illetve letiltásával a hagyományos NVT, NpT és µVT sokaságok mellett egyéb, tetszés szerint definiált sokaságokon is képes működni. A forráskódnak jelentős részét teszik ki a speciális szimulációs technikákat megvalósító rutinok (pl. az energia szerint irányított identitáscsere-eljárás [24] vagy a DCV-GCMD technikán alapuló módszer [28, 29]).

A zeolit adszorpciós vizsgálataival kapcsolatos eredmények

2. Elsőként metanol/víz és etanol/víz elegyek adszorpciós egyensúlyait határoztam meg NaA-4 zeoliton, nagykanonikus sokaságon végzett Monte Carlo szimulációkban, elsősorban az általam kifejlesztett NaA-4-potenciálmodell tesztelésének céljából [81].

3. A potenciálmodellek 378 K hőmérsékleten történt tiszta anyagokra és elegyekre vonatkozó részletes összehasonlítása azt mutatta, hogy míg a korábbi számításokban használt A modell kis és közepes nyomástartományokban a metanol vagy etanol adszorbeált mennyiségére jósol nagyobb értékeket a vízzel szemben, addig az új modell (B modell) - a kísérletekkel nagyobb összhangban - a vizsgált nyomástartomány egészében a vízre jósol nagyobb mértékű adszorpciót.

4. Az adszorbens adszorbátumhoz való affinitásának jellemzésére szolgáló zérus betöltöttségre extrapolált izosztér adszorpciós hő számított eredményei azt mutatták, hogy az A modellel megfigyelhető tendencia, miszerint az adszorbens-adszorbátum kölcsönhatások erőssége a víz, metanol, etanol sorrendben nő, az új modell használata esetén megfordul. A megfigyelések alapján arra következtethetünk, hogy az új modell viselkedése valószínűleg jobban megfelel a permeációs (nemegyensúlyi) kísérletekben mérhető nagy szelektivitásértékeknek, amelyek azt mutatják, hogy a NaA-4 zeolit lényegesen szelektívebb vízre a metanollal vagy az etanollal szemben.

5. A potenciálmodellek szerkezeti és energiaanalízise szintén alátámasztotta az új modell reális viselkedését. A szimulációkban meghatározott párkorrelációs függvények tulajdonságai erősebb adszorbens-adszorbátum kölcsönhatásokra utalnak a polárisabb új modell esetén.

6. Olyan sűrű rendszerek esetén, mint az NaA-4 zeolit, az adszorpció szimulációs vizsgálatát némelykor az adszorbátum kémiai potenciáljának meghatározása megkönnyíti.

Ehhez sikeresen használtam egy általunk fejlesztett, adaptív mintavételezéses termodinamikai integráláson alapuló módszert [87].

A kaolinit interkalációs vizsgálatával kapcsolatos eredmények

7. Munkám második részében a réteges szerkezetű kaolinit agyagásvány interkalációja során létrejött stabil komplexek jellemzésére használt d001 rétegtávolságokat, betöltöttségeket, illetve az interkalációs komplexek szerkezetét leíró eloszlásokat határoztam meg µVT, µpT és NpT sokaságokon végzett szimulációkban.

8. Megállapítottam, hogy minden egyes vizsgált interkalálószer (klasszikus interkalálószerek: DMSO, formamid, karbamid, kálium-acetát, továbbá a kísérletekben preinterkalátumok segítségével interkalálható anyagok: víz és metanol) esetén az interkalációs komplexeknek több, különböző rétegtávolsággal jellemzett stabil állapota detektálható a szimulációkban [124, 126, 128]. A jelenség vizsgálatára szimulációkat

végeztem az irodalomban kaolinitre [117, 118], illetve az interkalált molekulákra vonatkozó egyéb potenciálmodellekkel is (pl. [119]). A kapott eredmények alapján elmondható, hogy a két különböző szerkezettel és rétegtávolsággal jellemezhető stabil komplex létezése általános érvényűnek tűnik, és független az alkalmazott potenciálmodellektől. Esetleges további, nagyobb rétegtávolságú stabil komplexek létezésére a számításokban nincs bizonyíték.

9. A szimulációk során a különböző rétegtávolságoknál meghatározott szabadentalpia vagy rétegnyomás értékeinek változása azt mutatta, hogy többnyire egy meglehetősen nagy energiagát van a szimulációkban detektálható első és második stabil rétegtávolsággal jellemzett állapot között. Ez lehet az oka, hogy eddig kísérletileg nem sikerült igazolni azt, hogy a legalább két különböző rétegtávolságú stabil interkalációs komplex létezése általános érvényű lenne, noha a víz mellett kálium-acetátra vonatkozóan is közöltek olyan kísérleti eredményeket, amelyek valószínűsítik ezt a jelenséget.

10. A kálium-acetát/kaolinit interkalátumra vonatkozó szimulációs eredményeink szintén két termodinamikailag stabil állapotot mutattak, azonban az eddigi tapasztalatoktól eltérően a hagyományos interkalációs kísérletek a második stabil rétegtávolságú komplexet erősítették meg. A szimulációs eredményeink támpontot szolgáltattak egy új kísérletsorozat elvégzéséhez, amelyben megfigyelték a kisebb rétegtávolságú komplexet, és az egymást kölcsönösen támogató szimulációs és műszeres vizsgálatok tisztázták annak önálló, megfelelő körülmények között stabil voltát [128].

11. A szimulációk során meghatározott, a molekulák egyedi orientációját és a molekularétegek elhelyezkedését leíró eloszlásprofilok alapján megállapítható, hogy az interkalált molekuláknak az első stabil rétegtávolság esetén egy egyrétegű, míg a második (nagyobb) rétegtávolságnál egy kétrétegű elrendeződése alakul ki a kaolinitrétegek között, ahol az egyes molekularétegeken belül a részecskék többnyire egy jól meghatározott orientációt vesznek fel. Míg az egyrétegű elrendeződés esetén a molekulák síkja általában többé-kevésbé párhuzamos a kaolinitrétegekkel (bár a dipólusok lehetőség szerint enyhén az Si-O lemez felé mutatnak), addig a kétrétegű elrendeződés esetén az átlagában

párhuzamos elrendeződés csak diffúzabb módon érvényesül, sőt, a karbamidnál felfedezhetők közel merőleges molekulaorientációk is.

Dinamikus Monte Carlo szimulációk többkomponensű rendszerekre

12. Munkám harmadik részében időfüggő tulajdonságok sztochasztikus számítási lehetőségével foglalkoztam. Mivel a DMC szimulációban eltelt idő és ennek megfelelően a számított dinamikus tulajdonságok nem függetlenek a különböző szimulációs paraméterek (pl. rmax) értékeitől, ezért egy olyan eljárást dolgoztam ki, amellyel többkomponensű rendszerek komponenseinek relatív dinamikája helyesen számítható molekuláris dinamikai referenciaszámítások nélkül is, ha az adott rendszerben a transzláció a meghatározó mozgásforma [135]. A módszer alapvetően három fő kritérium köré csoportosítható:

• a hagyományos MC módszerhez hasonlóan a transzlációs szimulációs lépés részecskekiválasztási valószínűsége legyen 1/N;

• a fizikailag valószínűtlen lépések elkerülésére rmax értéke abból az átlagos úthosszból legyen meghatározva, amelyet egy molekula a szomszédos molekulákkal való ütközésig megtehet;

• a szimulációban kapott dinamikus tulajdonságokat skálázzuk át a komponensek tömegének négyzetgyökével.

13. A javasolt módszert több különböző rendszeren, köztük a membrántranszport direkt szimulációját megvalósító DCV-GCMD technikán alapuló számításokban teszteltem, összehasonlítva MD eredményekkel. Megállapítottam, hogy a javasolt módszerrel egy-egy rendszerre számított dinamikus tulajdonságok arányai az összes tesztrendszer esetén jó egyezést mutatnak az MD referenciaszimulációk megfelelő eredményeivel.

14. A DMC szimulációkkal vizsgált kalciumcsatorna (biológiai ioncsatorna) esetén kapott eredmények [145] összhangban vannak az eddigi kísérleti megfigyelésekkel abban a tekintetben, hogy a vizsgált ioncsatorna a teljes móltörttartomány túlnyomó részében a kalcium-ionra mutat nagyobb szelektivitást a nátrium-ionnal szemben.

15. A DMC szimulációs eredmények alapján elmondható, hogy ugyanolyan körülmények között a vizsgált ioncsatorna affinitását jellemző, betöltöttségek hányadosaként kifejezett szelektivitás a kalcium-ionra nagyobb, mint az ionáramok hányadosaként definiálható szelektivitás. Az általában nagyobb áramok ellenére ez azt jelenti, hogy a kalcium-ion nem szívesen hagyja el a csatornát. A bruttó betöltöttség és kalcium-ionáramok mellett számított molekulasebesség-profilok és betöltöttség-eloszlások elemzése alapján új megvilágításba kerülhet a kalciumcsatorna szelektivitása mögött rejlő mechanizmus.

SUMMARY

The most important results presented in this dissertation can be summarized in the following theses.

Thesis points:

Development of the simulation software

1. My initial aim was to expand the capabilities of the C program I previously developed.

With the platform independent, object-oriented new C++ version user defined ensembles can be realized, in addition to the traditional NVT, NpT and µVT ensembles, by enabling or disabling different types of MC moves. A large part of the source code contains routines that can handle special simulation techniques (e.g. the energy biased identity exchange technique [24] or the DCV-GCMD based method [28, 29]).

Results for the adsorption in zeolite

2. Equilibrium adsorption properties of water/alcohol mixtures in zeolite NaA-4 were determined by grand canonical MC simulations, primarily to investigate the behaviour of our novel (reparameterized) potential model for the zeolite [81].

3. A detailed comparison of the new potential model (model B) and the one I used earlier (model A) was made for pure components and mixtures by performing simulations at 378 K. This comparison showed that while model A predicts higher loadings for methanol and ethanol than for water in the low and medium pressure range, the new model predicts, in better harmony with experiments, higher loadings for water in the whole pressure range.

4. According to the simulation results, the strength of the adsorbent-adsorbate interaction, which can be characterized by the isosteric heat of adsorption extrapolated to zero

for the new model a reverse tendency can be observed. These findings suggest a more realistic behaviour of the new model considering that in permeation (non-equilibrium) experiments zeolite NaA-4 shows very high selectivity of water to alcohols.

5. The energy and structural analysis results also support the realistic behaviour of the new model. The properties of the calculated pair correlation functions suggest that the adsorbent-adsorbate interactions are stronger for the new, more polar model.

6. Determining the chemical potential for the adsorbate species can facilitate, in certain cases, the molecular simulation of the adsorption phenomena in dense systems. For the adsorption with zeolite NaA-4, I successfully applied our chemical potential calculation method that uses thermodynamic integration with adaptive sampling [87].

Results for the intercalation in kaolinite

7. In the next phase of my work, the characteristic properties (d₀₀₁ basal spacing and loading) of stable kaolinite intercalation complexes were determined in µVT, µpT and NpT ensemble simulations. The arrangements of the intercalated molecules were described by density profiles and orientation distributions.

8. In simulations with known intercalation compounds (DMSO, formamide, urea, potassium-acetate) and other substances that can be intercalated utilizing expansion precursors (water and methanol), two stable kaolinite intercalate complexes with different basal spacings can be detected [124, 126, 128]. Additional simulations were also performed with other potential models for kaolinit [117, 118] and for the intercalated molecules (e.g.

[119]) available in the literature. The results showed that the existence of two stable states seems general and does not depend on the applied potential models. There was no simulation evidence that supported the existence of additional stable states with larger basal spacings.

9. The tendency of the calculated Gibbs free energy and disjoining pressure calculated at different basal spacings suggests the existence of a large energy barrier between the two observed stable states, and this might be the reason that, generally, the second stable state has not been confirmed experimentally yet (although for water and potassium acetate there are some experimental results in the literature that seem to support this phenomenon).

10. In the case of potassium acetate, again, two stable states were detected in simulations, but in contrast to the other substances, the stable complex characterized by the larger d001

value was verified in experiments. Following the prediction of our simulations, new experiments and measurements were carried out, and the existence of the stable intercalate complex with smaller d₀₀₁ value was confirmed by a joint theoretical/experimental analysis [128].

11. The calculated density profiles indicated that two different interlayer structures are realized for the two detected stable states: depending on the basal spacing, either a single or a double layered arrangement of the intercalated molecules develops between the kaolinite layers. According to the orientation distributions the intercalate molecules have a well defined orientation between the kaolinite layers. For the single layered arrangement the molecular planes are approximately parallel with the clay layers (though dipole moment vectors tend to point towards the Si-O sheet), while for the double layered arrangement the orientation distributions become more diffuse and the parallel orientation is realized less frequently (for urea a more or less perpendicular orientation can also be identified).

Dynamic Monte Carlo (DMC) in many-component systems

12. The calculation of time dependent properties with the stochastic DMC method was the final task in my work. In DMC the time (and, therefore, dynamic properties) depends on the values of the simulation parameters (e.g. rmax). A new method was presented that allows the DMC procedure to realize the correct time proportionality without corresponding MD calculations in many-component systems, where translation is the dominant type of motion [135]. The method consists of three inherent parts:

• the standard MC selection probability 1/N for particles is used;

• the value of r_max is determined from the average path that a molecule can move towards its neighbours until collision to avoid physically unrealistic jumps in the DMC simulation;

• the calculated dynamic properties of every component are rescaled with the square root of the respective particle mass.

13. Comparing MC and MD simulation results, the method was tested in different systems including a test system in which a DCV-GCMD based method was applied that provides a possible solution to simulate membrane transport. The proposed method was found to provide consistently good results for the dynamic property ratios in each system as compared to the reference MD calculations.

14. The obtained DMC results [145] for a calcium channel (biological ion channel) are in agreement with the prior experimental observations: the channel shows higher selectivity for the calcium ion as compared to the sodium ion in the major part of the mole fraction

14. The obtained DMC results [145] for a calcium channel (biological ion channel) are in agreement with the prior experimental observations: the channel shows higher selectivity for the calcium ion as compared to the sodium ion in the major part of the mole fraction

In document Nanopórusos rendszerek vizsgálata molekuláris szimulációkkal (Pldal 97-123)