Molecular Dynamics Simulations to Study Drug Delivery Systems

Abstract

Molecular dynamics simulation is a very powerful tool to understand biomolecular processes. In this chapter, we go over different applications of this methodology to drug delivery systems (DDS) carried out in the group. DDS—a formulation or a device that enables the introduction of a therapeutic substance in the body and improves its efficacy and safety by controlling the rate, time, and place of release of drugs—are an important component of drug development and therapeutics. Biocompatible nanoparticles are materials in the nanoscale that emerged as important players, improving efficacy of approved drugs, for example. The molecular understanding of the encapsulation process could be very helpful to guide the nanocarrier for a specific system. Here we discuss different applications of drug delivery carriers, such as liposomes, polymeric micelles, and polymersomes using atomistic and coarse grain (CG) molecular dynamics simulations.

1. Introduction

Nowadays nanomedicines have a significant impact in healthcare, either for diagnosis or for therapeutic purposes, as in the case of drug delivery systems (DDS). These systems have been shown to protect the drug from external factors and/or degradation, to prolong the release rate, to target specific tissues, and to change organoleptic aspect [1, 2, 3]. The field of pharmaceutical research has developed enormously in recent years, providing new DDS for the transport and release of bioactive compounds. To achieve that, novel nanomaterials have improved, thanks to the rapid advances in material sciences, many of them with tunable properties that allow the controlled release of drugs [4]. Such nanoparticulated systems include polymeric micelles, polymer-DNA complexes, nanogels, lipid-based (liposomes, nanostructured lipid carriers) substances, and macrocyclic carriers (cyclodextrins, calixarenes), among others. DDS provide an extraordinary opportunity for the safe and efficient release of drugs, genes, and a great variety of molecules. The incorporation of nanomaterials confers physical advantages such as improved drug solubility, decreased degradation or clearance rates, decreased systemic toxicity, and improved clinical efficacy to the nanoformulations [5]. Furthermore, the sustained release of the drugs allows one to reduce the frequency of drug administration or even the dose. Finally, the process time for approval of a novel DDS that carries a drug already sanctioned by the regulatory agents is shorter, making the development process less expensive.

Relevant information on the system can be obtained using computer simulations at very different lengths and time scales: moving from continuous to the atomistic level. In this context, knowledge of the mechanism of drug encapsulation and release at the atomic/molecular level can help in the design of nanomedicines, according to the desired objective and for each specific case.

In this sense, techniques that access the molecular level, such as molecular dynamics (MD) and Monte Carlo simulations, are very powerful tools to understand biomolecular processes [6, 7]. In particular, classical MD simulations could help in the development/improvement of drug delivery systems. But, at the current computer capacity, not all nanoparticles can be studied by MD, at their full length. How to approach to this problem? For instance, liposomes are widely used DDS, consisting of spherical lipid-based vesicles whose diameters range from 30 nm to several micrometers. To simulate this kind of DDS carrier there are two possible approaches, depending on the question to answer: (i) considering a big liposome, that is, having an infinite radius, one can use planar bilayers or also simulate just a section of the vesicle, within periodic boundary conditions and (ii) on the other hand, a small liposome can be fully treated. This is exemplified in Figure 1. Broadly, we can access atomic details of drug-bilayer interaction using atomic level MD, but liposomes are better signified using coarse grain (CG) models.

Figure 1.

Schematic representation of the section of the vesicle bilayer and small liposome used as model for atomistic and coarse grain simulations.

In the following sections we will discuss fundamentals of the MD simulations and examples of simulations of drug delivery systems at the atomistic and CG level. A brief description of some approaches to study the release of drugs is also included. We finish the chapter with, in our opinion, important open questions in the field.

2. Methodology

The molecular dynamics technique has a long history [8] and it nowadays constitutes a very important theoretical tool in physics, chemistry, biology, and related disciplines. Thanks to the growing development of computing power, the MD technique is renewed continuously, allowing one to study larger and more complex systems and problems [8].

2.1. Classical mechanics approach

The physical properties of matter are associated with its structure and the movement of its basic constituents, nuclei and electrons. In principle, the evolution of a system of particles—a many-body problem—can be obtained by solving the Schrödinger equation dependent on time, which gives the probability of finding the particles in any position of space in a given time. In many cases, it is possible to decouple the wave function of the electrons from the wave function of the nuclei, following the Born-Oppenheimer approach (due to mass difference). However, even considering this approach, it is impossible, in practice, to solve the Schrödinger equation numerically and further approaches must be used (Hartree Fock, density functional theory, etc.). In particular, for biological systems (where thousands of nuclei and electrons), its resolution is impossible even under such approximations.

In this direction, molecular dynamics simulation could help to study problems of many bodies at the atomistic level, based on classical mechanics. Within this approach, Newton's equations are solved numerically. The main advantage is the realistic simulation of materials through the simplification by potentials with analytical form.

Although it is quantum mechanics, instead of classical mechanics, which describes the fundamental physics of condensed matter, the validity of the classical approximation can be evaluated based on the Broglie thermal wavelength defined by

$$\mathrm{\Lambda }={\left(\frac{2\pi {\mathrm{\hslash }}^2}{{\mathit{mk}}_{b}T}\right)}^{1/2}$$E1

where ħ is the Plank constant, m the atomic mass, k_b the Boltzmann constant, and T the temperature. The classic approach is valid for Λ < <a, being the separation of the nearest neighbor. Under this condition, the entire system can be treated as a diluted gas model based on the formulation of classical kinetic theory [9]. In this case, each atom can be considered as a particle. Many biomolecular processes can be addressed by semi-empirical parameterizations that describe the interactions of pairs between the particles of the system with classic effective potentials. In addition to having a lower computational cost, the classic models offer an adequate description of the processes and interactions. Chemical reactions (formation and breaking of bonds) and vibrations of bonds and angles with very high frequencies are excluded from the classical treatment.

Through MD simulations, information of positions and velocities of each particle is obtained, from which macroscopic observables could be calculated, such as pressure and energy.

2.2. Statistical mechanics: averages in a simulation

Here, we would like to discuss the relationship between observed properties of a large system and their microscopic dynamics or fluctuations. The particles of interest (atoms, molecules, or electrons and nuclei…) obey certain microscopic laws with specific interparticle interactions. The task of solving the equation of motion for a many-body system is still complicated and difficult, even with nowadays computer power.

Experiments are usually done on a macroscopic sample, containing a number extremely large of atoms or molecules in a huge number of conformations. In statistical mechanics, the average of the observable A over a given ensemble is calculated on a large number of replicas of the system to obtain the observables of the experiment.

On the other hand, the microscopic state of a system is defined by the positions and velocities of the particles, which are the coordinates of a 6 N multidimensional space (where N is the number of particles). The experimental observables are the averages of the ensemble and not the temporal averages. However, in MD simulation, thousands of atoms/beads are used as a way of sampling of a mechanical-statistical ensemble. How to solve this difference between the temporal averages and ensemble averages? The answer leads us to one of the fundamental axioms of statistical mechanics, the ergodic hypothesis, which establishes:

$$⟨\mathrm{A}⟩\text{ensemble}=⟨\mathrm{A}⟩\text{time}.$$E2

By allowing the system to evolve indefinitely, this means that it would pass through all possible states compatible with the constraints. Even if it impossible in practice, it is important to guarantee a wide sampling of representative conformations of the phase space when MD is used.

The ensembles are characterized by constant values of thermodynamic variables that describe the state of the system. Thinking of the simulations as a computational experiment, the statistical ensemble characterizes the conditions in which an experiment is performed. For instance, NPT ensemble has constant temperature (T), pressure (P), and number of particles (N). Each particular state, defined by these parameters, has an associated state equation that characterizes the system.

In order to mimic biologically relevant model membranes, many simulations have been done considering positive surface tension (ensemble NP_zγT). This ensemble is chosen due to the fact that experimental bilayers are able to adjust their area per lipid in order to achieve a minimum of free energy. However, it is argued that in simulations, periodic boundary conditions limit the bilayer undulations that impact the interpretation of the surface tension. In this way, full relaxation of the simulation box (NPT ensemble) is the most used for these kinds of systems [7, 10, 11, 12, 13].

2.3. Thermostat and barostat

In a molecular dynamics simulation, integrating Newton's equations of motion provides the means to sample the physical characteristics of a given system through its evolution in the microcanonical ensemble (NVE) [14]. Nevertheless, in order to sample other ensembles, additional variables should be added, known as extended degrees of freedom. For instance, to keep a constant temperature, the additional variables are used to control the temperature through the use of a “thermostat” [15, 16, 17]. In a similar way, in order to maintain a constant pressure during the simulation, a “barostat” constructed with additional pressure-controlling variables should be added [18, 19]. In this way, the average temperature and pressure are regulated within these variables and their dynamics. Aiming to achieve such temperature and pressure regulation without disturbing the short-term, Newtonian dynamics of the particles, the extended degrees of freedom are generally designed to evolve slowly and with only weak coupling to the dynamics of the physical particles. Whereas the fastest timescales of atomic motion are in the tens of femtoseconds, the timescales for the extended degrees are typically picoseconds or longer [14].

2.4. Boundary conditions

To perform an MD simulation, it is necessary to define the “simulation box” containing different molecules, subject to the appropriate boundary conditions to the geometry of the macroscopic system. The number of particles in the simulation box is very small, compared to an experimental sample. In this context, the border effects are meaningful. In order to minimize, usually, periodic boundary conditions could be considered. In this approach, the simulation box is surrounded by identical copies in all directions, giving rise to a periodic system that tends to the thermodynamic limit. When a particle leaves the simulation box, its image enters simultaneously through the opposite face. For bilayers, this kind of conditions is necessary in order to carry out the simulations. For other systems, like liquid droplets, non-periodic boundary conditions could be used [20].

2.5. Level of description and force fields

As we already mentioned, here we will discuss two description scales within a classical approach for the systems: atomistic and coarse grain. The main difference between them is that in an atomistic scale all the atoms (even hydrogen) are represented, whereas in a coarse grain system, atoms are grouped in beads [21]. In this way, CG models allow not only the reduction of the degree of freedom but the possibility to integrate Newton equations in a higher time step, due to the elimination of high-frequency vibration modes [22]. Usually, the treatment of interactions in the simulated systems involves the introduction of an effective force field, which allows performance of large-scale calculations of relatively large membranes and nanoparticles. This set of interaction models contains all functions and values of their parameters for the simulation.

In force fields, the potential energy could be divided into two parts, the intramolecular potential energy (the bonded part) and the intermolecular potential energy (the non-bonded part). Bonded atoms usually involve the stretching, the bending, and the torsion of the bonds. The non-bonded potential energy could be modeled as Lennard-Jones interactions and electrostatic interactions. These models are also applied between atoms of the same molecule, which are more than three bonds apart. Atoms that are three bonds apart (1–4 pairs) are treated differently in each force field.

The force field parameters are fitted to experimental and quantum mechanical data to match the spectroscopic, thermodynamic and crystallographic data of the molecule or, in the case of CG force field, they could also be derived from atomistic simulations. Among the classic effective potentials available are the atomistic force fields OPLS [23], CHARMM [24], AMBER [25], GROMOS [26] and CG Martini [21], SIRAH [27], and others [28]. The choice of the force field is a crucial step, and it should be carefully chosen since it strongly affects the results of an MD, and the number of works are constantly checking the precision of them [29].

3. What we can get from simulations?

3.1. Atomistic simulations

Within the fully atomistic treatment, only few nanocarriers could be completely simulated. This is the case, for instance, for cyclodextrins [30, 31], calixarenes [32], dendrimers of low numbers [33], and so on. Nevertheless, in most cases, in order to get information of the nanoparticles, further approaches should be done:

• small models of the complete systems: small micelles, nanostructured particles [34];

• focus on single-drug interaction with specific components of the system [35];

• internal structure of the nanoparticle in the presence of the drug;

• bilayers as a model of liposomes, polymersomes, niosomes, etc. [36, 37].

Is in this last approach where we will emphasis in this section due to the amount of work present in the literature [7, 38, 39, 40, 41, 42, 43] and the direct connection with experimental model bilayers. Within atomistic simulations it is possible to explore, for instance, the drug partition in the different bilayer regions and go deeper into the main interactions responsible for their localization [37]. The drug localization could be guided by specific interactions, such as hydrogen bond [36], cation-π (for aromatic molecules) [44], salt bridges [45, 46], or entropic effects (hydrophobic drugs usually partitioned in the hydrophobic core [47].

Several biophysical techniques could be used to access bilayer properties. Between them, X-ray and neutron-scattering experiments allow us to obtain the electron density profiles that give us an idea on the overall organization of a lipid bilayer. Within the simulation not only it is possible to reproduce this profile but also access each system component contribution. The profiles are calculated dividing the simulation box in slices normal to the bilayer and conducting the time averaging of the net charge (equivalent to the total electron associated with each atom) per given thick slabs. Other quantities like radial distribution function could be used in order to obtain more details on specific interactions between the drug and the lipids or to extract information about water networks and compare with the structure factor of X-ray experiments [20].

Local anesthetics (LA), pain-relief drugs, are among the most extensive drugs studied by simulations in their interaction with membranes. The structure and physicochemical features of each LA agent determine its potency and toxicity. Encapsulation into DDS -such as liposomes- can modulate the anesthetic effect (onset of action, duration of sensory block) and toxicity [5]. In this direction, simulations could guide the development of DDS devoted not only for short-term procedures (ambulatory and surgical interventions) but also for chronic pain treatment.

Most LAs have a pKa between 7 and 9; therefore, both neutral and protonated forms are relevant at physiological pH. Most works in literature overcome this issue studying the drug at a high and low pH, where all molecules are protonated or deprotonated, respectively. From MD simulation, most authors report neutral LA partition in the hydrophobic region of the bilayer and the protonated species in the lipid head/water interface and water phase, with differences between them (depending on hydrophilicity, steric factors, the lipid bilayer, etc.) [47, 48, 49, 50, 51, 52, 53]. For instance, the behavior described for lidocaine in 1,2-dimyristo yl-sn-glycero-3-phosphocholine (DMPC) is similar to the one reported by us, using prilocaine in 1-Palmitoyl-2-oleoylphosphatidylcholine (POPC), even using different force fields (GROMOS and CHARMM27, respectively) [48, 50]. In both cases the neutral drug presents a bimodal distribution. Also, crossing events, defined as the drug crossing from one monolayer to the other bilayer, were observed for both drugs. On the other hand, in our knowledge, no crossing events were reported for any protonated LA, even in long simulations (400ns, Slipid force field ) [53]. The differential distribution of prilocaine (PLC) depending on its ionization state is illustrated inside the bilayer as shown in Figure 2A. On the other hand, more hydrophobic LA at the neutral state, such as etidocaine, shows higher distribution than prilocaine in the bilayer core [50]. The localization of the etidocaine in this region promotes a higher disorder in the lipid chains. This property could be accessed through the calculation of the order parameter, defined by:

$$S_{\mathit{mol}}=\frac{1}{2}⟨3{\text{cos}}^2{\theta }_n-1⟩$$E3

Figure 2.

(A) Electron density profiles (EDPs) of different bilayer components as a function of the membrane normal: POPC in black, water in blue and drug in red. Solid and dotted line are used for protonated and neutral prilocaine, respectively. (B) Calculated order parameter, −SCD, for plain bilayers (black); etidocaine at 1:6 and 1:3 drug:lipid ratio in green and violet, respectively. Circles and circles are used for saturated and unsaturated chains of POPC.

Here, θ_n is the angle between the normal to the bilayer and the normal to the plane defined by two carbon-hydrogen bonds. The order parameter is related to the tilt angle of the chains and to the trans-gauche distribution of chain dihedrals. The effects on the lipid packing from guest molecules such as PLC distributed within membranes can be experimentally determined from changes in the carbon-deuterium segmental order parameter along the lipid chain. In that case, the experimental order parameter, S_CD = -1/2S_mol, is derived from the measured quadrupole splitting Δν = (3/4)(e²qQ/h)S_CD [54, 55]. In Figure 2B we show the order parameter of POPC neat and containing etidocaine.

In a similar direction, simulations with the two bupivacaine (BVC) enantiomers helped to elucidate the higher cardiotoxicity observed for the R-form than the S-form [56]. The effect of R-BVCs (at 1:3 molar concentration) is to disorganize the membrane (decrease the order parameters); this effect is seen for both tails [51]. This is essentially related with the empty space in the lower part of the lipid tails due to the localization of the LA and the lateral expansion of the bilayer. On the other hand, S-BVC only promotes a soft disorganization for carbons above 10 (in both tails): in this case, the localization of the LA in the interior of the bilayer compensates the effects of the lateral expansion.

Thus, in our hands, molecular simulation of LA in model membranes was found useful to explain different aspects on the anesthesia mechanism and drug encapsulation [57]. Besides, within simulations, it was possible to explain the experimentally observed differences [49] between less (prilocaine/lidocaine) and more hydrophobic (etidocaine) isomers regarding the depth of their preferential insertion into bilayers [50], with possible implication on the increased potency and toxicity of the more hydrophobic analogs (Meyer-Overton rule) [5, 58].

3.2. Coarse grain

Through the reduction of degrees of freedom, CG models are useful to efficiently simulate drug delivery systems, such as liposomes, polymersome, and micelles, relieving the size and time-scale limitations of atomistic simulation but losing in details. One of the crucial factors, representing the capacity of potential drug delivery systems, is the partition coefficient of a potential drug candidate between the aggregate and surrounding water.

In particular, the encapsulation of prilocaine into liposome was studied using this approach (MARTINI-based force field) [59]. Following the atomistic results, neutral PLC was fully encapsulated in the interior part of the lipid membrane where it adopts an asymmetric bimodal distribution. Our simulation results therefore suggest that although protonation leads to a structured interaction between drug and host, hydrophobicity is the major driving force of drug encapsulation. Our results also depend on the protonated PLC initial simulations conditions [59]. This observation suggested different preparation schemes of liposome-drug complexes, leading to prilocaine trapped within vesicles that could increase overall drug encapsulation efficiency.

On the other hand, it is demonstrated that the two LA species (neutral and protonated) are present at physiologic pH, contributing together to the anesthetic effect [50, 57] and they could be important in the development of DDS. For example, in a liposomal formulation prepared at pH 7.4, while the neutral form of the anesthetic is able to quickly and efficiently cross the bilayer, the protonated form is mainly found in the water phase. Both species depend on each other: the protonated needs the neutral one in order to efficiently cross the membrane, reaching the adjacent water compartment, while the uncharged species depends on the high solubility of the protonated LA in the water phase to reach the clinically effective dose.

The protonation/deprotonation reaction cannot be directly simulated using classical MD simulation. In order to partially overcome this, we represented the physiological pH, taking into account Henderson-Hasselbach equation:

$$\mathit{pH}=\mathit{pKa}+\text{log}\frac{\left[\mathit{\text{neutral}}\mathit{PLC}\right]}{\left[\mathit{\text{protonated}}\mathit{PLC}\right]}$$E4

and the fact that PLC apparent pKa in membranes is 7.6 [60]; values around 0.4 and 0.6 represent the molar fractions of neutral and protonated species, respectively.

In Figure 3A, we illustrated the results with a snapshot. In order to get an idea of how PLCs distribute inside the vesicle, we have calculated the average number of PLCs as a function of the PLC (center of mass) distance to the vesicle center. In Figure 3, we also show separated histograms for neutral (B) and protonated (C) PLCs. Considering that the average vesicle radius is 75Å, we can see from Figure 3C that all neutral PLCs are essentially found inside the vesicle and show two main peaks, in good agreement with previous results [47, 57]. On the other hand, just a small fraction (~14 molecules on an average) of the protonated PLCs are found inside the vesicle, as we can see in Figure 3C, and they only interact with the external monolayer of the vesicle. The results here show that the behavior of PLCs at physiological pH is essentially a combination of high and low pH: This means neutral PLCs are found inside the vesicle, whereas protonated molecules are partitioned into the external monolayer of the vesicle and water regions. However, different than at low pH, the average number of protonated PLCs inside the vesicle is lower at physiological pH (14 than 20) because of the presence of neutral molecules. In this direction, this kind of simulation could guide the LA:lipid molar ratio to enhance the loading capacity of the liposomes, depending on pH.

Figure 3.

(A) and (B) Snapshots from the molecular dynamics simulations prilocaine in POPC vesicles at physiological pH. The snapshots correspond to the initial configuration and after 1 μs NVT simulation run. Neutral and protonated PLC molecules are shown in magenta-blue and in brown, respectively, and lipids are in blue and green forming the vesicle. Water sites were removed for visualization purposes. The chlorine counter ions, in red, are only shown in the last configuration. Number of neutral (C) and protonated (B) PLCs as function of the PLCs center of mass distance to the center of the vesicle.

Similar to liposomes, polymersomes form a core-shell structure allowing us to encapsulate both hydrophilic and hydrophobic molecules. These synthetic polymer vesicles are composed of block copolymer amphiphiles that in the presence of water self-assemble in an organized structures, attracting special interest due to their tunable properties for drug delivery. However, since polymersome sizes often range from 100 nm up to the micron diameter size it is difficult to simulate, even at the coarse grain level, a big polymersome. In this direction, the simulation of polymeric bilayers and small polymersomes at a CG level would make possible to study key mechanical and structural properties to better develop these kinds of drug delivery systems, as recently explored by us [12]. Moreover, these simulations are useful to shed light into the effects of the incorporation of model drugs (such as neutral and protonated prilocaine), as we explore in a work-in-progress study.

Other popular drug delivery systems composed of copolymers are micelles. Several examples are found in literature [61, 62, 63, 64, 65, 66]. In this regard, the encapsulation of the hydrophilic antimigraine drug sumatriptan in a polymer micelle was studied by us in a recent work [13]. A micelle composed of the pluronic F127 tri-block copolymer was simulated under different initial conditions. The main results showed that the drug essentially partitioned in the hydrophilic drug with little effect in the overall micellar structure, and size also correlated with experimental studies.

3.3. Drug release: free energy calculations

The times required to observe the mechanism of drug release are out of reach using MD simulations. Nevertheless, important information could be obtained by calculating, for instance, the free energy through a given path for the drug leaving the nanostructure. The idea of enhancement of the sampling of configurational space is not new in molecular simulations, and several methods such as umbrella sampling [67, 68], adaptive biasing force (ABF) method [69], the Wang-Landau algorithm [70], steer molecular dynamics [71], and metadynamics have been proposed [53]. In order to investigate the release of the drug using these methods, a reaction path should be chosen. For bilayers, the natural path is the z direction (normal to the bilayer) and for spherical nanoparticles the radius of it, as exemplified in Figure 4A and B, respectively.

Figure 4.

Schematic representations (A) a free energy profile over the lipid bilayer normal and (B) a micelle, as generic DDS, showing the reaction pathway.

For the case of local anesthetics, very recently, Saaedi et al. estimated the free energy profile of lidocaine and articaine in a DMPC lipid bilayer using well-tempered metadynamics simulations (Slipid force field) [53]. They estimated that the free energy between the well and the water phase was −32.9kj/mol and −25.4 kj/mol for neutral lidocaine and articaine, respectively. On the other hand, this value reduces ~20% for the protonated species. Similar results were obtained by Prates et al. (unpublished results) using ABF, who reported a free energy difference ~24kj/mol in a POPC membrane using the NP_zAT ensemble. This means that pH is a determinant factor for the encapsulation/release of these drugs.

These calculations could be also carried out using a CG approach. Nevertheless, qualitative information could come across of these kinds of calculations. For example, Loverde et al. used Steered molecular dynamics to estimate the free energy profile of Taxol (an anticancer drug) when pulled out from a micelle core [72].

4. Final remarks and challenges

The rapid development of the field due to the entanglement between computer capabilities and algorithm development, open the possibility to imagine many applications of the MD approach in drug delivery systems. Till now, we have explored different applications on the molecular dynamics techniques to study drug delivery systems based on the available literature.

We would like to mention what we consider to be the next desirable steps in the field:

• In our knowledge, most of free energy calculations reported in the literature were simulated with a single drug, without taking into account the drug-loaded medium. In this way, important information could be obtained from these calculations.

• To obtain semi-quantitative information into the loading capacity of a specific drug into a given DDS.

• The use of a constant pH ensemble [73] could be applied to these kinds of systems. This represents a challenge in order to get on-the-fly protonation/deprotonation of molecules depending on the surrounding environment. This will emulate the loading and release of the DDS that could happen in different conditions.

• Through the calculation of the different components of the pressure tensor, it’s possible to access mechanical properties that could guide the DDS development.

• New DDS systems could be addressed with this technique. One example of this is the reconstituted high-density lipoprotein particles to deliver hydrophobic drugs to impaired cells and tissues [74]. These nanostructures are being used for other things such as protein stabilization, but the possibility to develop drug delivery systems with them is a promising one that could take advantages of molecular dynamics simulations.

• New applications for old drugs, as well as new therapies like gene delivery, could benefit from this computational approach [75].