complex chemical dynamics through engineering-like - AMS Dottorato

complex chemical dynamics through engineering-like - AMS Dottorato

` di Bologna Alma Mater Studiorum · Universita DOTTORATO IN SCIENZE CHIMICHE Ciclo XXVI Settore Concorsuale di A↵erenza: 03/A2 Settore Scientifico Dis...

17MB Sizes 0 Downloads 4 Views

Recommend Documents

tesi dottorato defitivita - AMS Dottorato
GABRIELE BALDUCCI. Prof. FABIO PEZZI ..... del prodotto (dato dal fatto della grande quantità di prodotto ammostato) e

Schema frontespizio teso Dottorato - AMS Dottorato
List of Figures. 6 .... African Union‟s Inter-african Bureau of Animal Resources ..... forced by destitution to depend

chapter 1 - AMS Dottorato
satellite images for flood mapping and monitoring, damage assessment and risk assessment. The contribution of satellite

Tesi Dottorato Stanzani-Maserati + FRONTESPIZIO - AMS Dottorato

Schema frontespizio teso Dottorato - AMS Dottorato
Emma, Middlemarch, David Copperfield. 3.0. Introduzione p. 137. 3.1. Le traduzioni italiane di Emma (1816) p. 143. 3.1.1

Facsimile del frontespizio della tesi di dottorato - AMS Dottorato
Mar 11, 2015 - of oil is equivalent to ten minutes of operation at low power or three ...... sidering the average per ca

engineering agent-oriented technologies and - AMS Dottorato
II Engineering Agent-Oriented Technologies for Programming Multi-. Agent Systems .... 7.7.1 Comparison with State-of-the

Facsimile del frontespizio della tesi di dottorato - AMS Dottorato
List of ..... of Molecular Biology, 1974-1980” Isis 92.3 (2001): 541-75; Arnold Thackeray (ed.) ... the grocery store,

ingengeria elettronica, informatica e delle - AMS Dottorato
Systems-on-a-Chip (SoC), the design space exploration needed to find the best .... Contents i. List of Figures v. List o

Low-temperature thermochronological evolution of - AMS Dottorato
Low-temperature thermochronological evolution of the Menderes and Alanya massifs. (Turkey). Presentata da: Dott. Frances

` di Bologna Alma Mater Studiorum · Universita DOTTORATO IN SCIENZE CHIMICHE Ciclo XXVI Settore Concorsuale di A↵erenza: 03/A2 Settore Scientifico Disciplinare: CHIM/02


Tesi di Dottorato in Chimica Fisica Presentata da: Lorenzo Moro

Relatore: Prof. Francesco Zerbetto

Coordinatore Dottorato: Prof. Aldo Roda

Anno Accademico 2012-2013 Esame Finale anno 2014


A quei due matti dei miei genitori, Alla mia famiglia, Ai miei amici, a Lollo (e anche a Marco, inestimabile compagnia in montagna), A Gilda e alla Marmolada. E al professor Zerbetto, a cui voglio bene. E anche alla bicicletta.


Contents Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction

vi 1

An engineering approach: a brief history . . . . . . . . . . . . . . . . . . . . .


Modern Engineering . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Computational Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . .


Rationale of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


2 Partial Di↵erential Equations: applications and solution


Partial Di↵erential Equations; who and why . . . . . . . . . . . . . . . . . . .


Boundary conditions and Initial Conditions . . . . . . . . . . . . . . . .


Numerical methods for solving di↵erential equations . . . . . . . . . . . . . .


Finite Di↵erences method . . . . . . . . . . . . . . . . . . . . . . . . . .


Finite Element method . . . . . . . . . . . . . . . . . . . . . . . . . . . .


COMSOL Multiphysics FEM environment . . . . . . . . . . . . . . . . . . . .


References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


3 PDE to study the molecular motors: the Fokker-Planck equation


Natural and artificial molecular motors . . . . . . . . . . . . . . . . . . . . . .


Rotary molecular motors . . . . . . . . . . . . . . . . . . . . . . . . . . .


Under the hood: The physics of Molecular Motors dynamic . . . . . . . . . .


Brownian Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Biased Brownian Motion . . . . . . . . . . . . . . . . . . . . . . . . . . .


A stochastic approach: the Fokker-Planck equation . . . . . . . . . . . .


References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .



4 Rotary molecular motors dynamic revealed


Theoretical framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


The chosen molecular motor . . . . . . . . . . . . . . . . . . . . . . . . . . . .


PES construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Escape rates, friction coefficients and Solution of the Smoluchowski equation


Results and analysis of them . . . . . . . . . . . . . . . . . . . . . . . . . . .


Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


5 PDE and nervous signals: the Hodgkin and Huxley model


Physiology of an excitable cell: an overview . . . . . . . . . . . . . . . . . . .


Electrical activity of the cell membrane . . . . . . . . . . . . . . . . . .


Membrane electrical properties: Calculating the membrane potential . .


Membrane electrical properties: The action potential . . . . . . . . . . .


Membrane electrical properties: Action potential propagation . . . . . .


Membrane electrical properties: Experimental methods for measuring action potentials . . . . . . . . . . . . . . . . . . . . . . . . . . .


Voltage Clamp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Target: separation into individual ionic components . . . . . . . . . . .


Voltage Clamp measurements . . . . . . . . . . . . . . . . . . . . . . . .


NA/K ionic currents separation . . . . . . . . . . . . . . . . . . . . . . .


The Hodgkin and Huxley model . . . . . . . . . . . . . . . . . . . . . . . . .


Potassium conductance: . . . . . . . . . . . . . . . . . . . . . . . . . . .


Sodium conductance: . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Transfer rate coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . .


Extensions of the HH model . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Membrane action potential . . . . . . . . . . . . . . . . . . . . . . . . .


Stimulation of excitable tissues . . . . . . . . . . . . . . . . . . . . . . .


Propagating nerve pulse . . . . . . . . . . . . . . . . . . . . . . . . . . .


Inclusion of the temperature e↵ects . . . . . . . . . . . . . . . . . . . . .


References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .



6 FEM simulation of Action Potentials


Validation of our FEM environment . . . . . . . . . . . . . . . . . . . . . . .


Numerical simulations of the Voltage Clamp . . . . . . . . . . . . . . . .


Membrane Action Potential Generation . . . . . . . . . . . . . . . . . .


Stimulation of active fibres with current pulses. . . . . . . . . . . . . . .


Propagating nerve pulse . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Computing Extracellular action potentials . . . . . . . . . . . . . . . . . . . . 105 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Model coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Computation of the extracellular potentials. Model validation and results.108 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7 Discussion and Future Perspectives


References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119


Glossary PDE - Partial Di↵erential Equation. ODE - Ordinary Di↵erential Equation. FDM - Finite Di↵erences Methos; The simplest numerical approach to solve PDEs. FEM - Finite Element Method; Extreme powerful numerical method used for obtain solutions of PDEs which cannot be solved analitically. FEA - Finite Element Analysis; Referred to a particular software which uses the FEM as a method to simulate the behaviour of system described by PDEs. CFD - Computational Fluido-Dynamic; branch of fluid mechanics that uses numerical methods and algorithms to solve and analyze problems that involve fluid flows so described by the Navier-Stokes equations. FPE - Fokker-Planck Equation; A particular PDE which describes the spatio-temporal evolution of the probability density function of the velocity or position of a particle population under the influence of drag forces and random forces, as in Brownian motion. HH - Hodgkin and Huxley; in this manuscript is used as HH Model referring to the model which explain the ionic mechanisms underlying the initiation and propagation of action potentials in the squid giant axon described for the first time in 1952 by the two scientists Alan Lloyd Hodgkin and Andrew Huxley.



Introduction Most of the problems in modern structural design can be described with a set of equations; solutions of these mathematical models can lead the engineer and designer to get info during the design stage. The same holds true for physical-chemistry; this branch of chemistry uses mathematics and physics in order to explain real chemical phenomena. With a typical engineering approach is thus possible to study the mathematical model related to a chemical process in order to fully understand and extract a lot of data which can help scientists in the development of new, more efficient, chemical structures.

An engineering approach: a brief history In the middle ages and throughout the ancient history, structural engineering was an art rather than, as today, a science; most of the architectural design and constructions were carried out by artisans, such as stone masons and carpenters, rising to the role of master builder. No record exists of any rational consideration, either as to the strength of structural members or as to the behaviour of structural materials; there were no theories about structures and all the understanding of how they stood up was extremely limited and based entirely on the empirical concept of what have worked before. Knowledge and experience of the builders were passed from generation to generation, guarded by secrets of the guild, and seldom supplemented by new knowledge. The result of this is that structures were repetitive, and increase in scale were incremental [1].


The first phase of modern engineering emerged in the Scientific Revolution. Galileo’s Two New Sciences [2], which seeks systematic explanations and adopts a scientific approach to practical problems, is a landmark regarded by many engineer historians as the beginning of structural analysis, the mathematical representation and design of building structures. Galileo’s studies were followed in 1676 by Robert Hooke’s first statement of the Hooke’s Law, providing a scientific understanding of materials elasticity and their behaviour under the presence of loads [3]. These studies were followed in the 17th century by Sir Isaac Newton [4] and Gottfried Leibniz which independently developed the fundamental theorem of calculus, providing one of the most important mathematical tool in the new born science engineering. From these milestones in the development of mathematical models which can helpful the design and construction of real world structures, we arrive to the modern era in which almost everything can be designed, tested and inspected before the real construction of it; those are the times of the modern engineering and of the computer aided technologies.

Modern Engineering At the present times, engineering mean the application of scientific, economic, social and practical knowledge in order to successfully design, build, maintain and improve structures, machines, devices, systems, materials and processes. In the last four thousands years of human evolution, engineering becomes extremely broad, and encompasses a range of more specialized fields, each with a more specific emphasis on particular areas of technology and types of application. The American Engineers’ Council for Professional Development has defined ”engineering” as: ”The creative application of scientific principles to design or develop structures, machines, apparatus, or manufacturing processes, or works utilizing them singly or in combination; or to construct or operate the same with full cognizance of their design; or to forecast their behaviour under specific operating conditions; all as respects an intended function, economics of operation or safety to life and property.” [5]


The basic engineer methodology lies on the creation of an appropriate mathematical model to find suitable solutions to a problem in order to analyse it and to test potential solutions. As with all modern scientific and technological endeavours, computers and software play an increasingly important role; there are a lot of computer aided applications (knowed as computer-aided technologies) developed specifically for engineering purposes. Computers can be used to generate models of fundamental physical processes, which can be solved using numerical methods. An example of how this tools can be successfully applied is represented in figure 1.1:

Figure 1.1: - Computer simulations of heat flow around the space shuttle during the re-enter in earth atmosphere. Credits NASA.

In the upper figure a computer simulation of high velocity air flow around the space shuttle during the atmosphere re-entering is showed. Solutions to this flow require modelling of the combined e↵ects of both fluid flows and the heat equation. With the modern computers is thus possible to fully simulate the behaviour of a complex structure or machinery so that is possible to fully design, test and inspect a construction before to really construct it. One great example of this approach is represented by the Virgin Formula one race car which have been completely designed, tested and developed in a virtual environment [6].


Computational Chemistry The science which applies mathematics and physics to chemistry with the aids of computers is the computational chemistry which is a branch of chemistry that uses computer simulation to assist in solving chemical problems. It uses methods of theoretical chemistry, incorporated into efficient computer programs, to calculate the structures and properties of molecules and solids [7]. It is widely used in the design of new drugs and materials. As well as for human dimension problems, many fundamental laws of physics and chemistry can be formulated as di↵erential equations. As an example many di↵usion processes can be described by second order partial di↵erential equations from whose solutions one can extract extremely useful informations. Studying chemical processes by this approach can complement the information obtained by real experiments and it can, in some cases, predict unobserved chemical phenomena. In this work we want to go beyond the traditional computational chemistry tools and apply a typical engineeringlike approach in order to describe and fully understand two extremely di↵erent chemical and physical processes.

Rationale of this thesis In this thesis work we will study two extremely di↵erent chemical processes: firstly we will study the dynamic of an artificial molecular motor and then the generation and propagation of the nervous signals between excitable cells and tissues like neurons and axons. These two processes, in spite of their chemical and physical di↵erences, can be both described successfully by partial di↵erential equations, that are, respectively by the Fokker-Planck equation [8] for the dynamic of the molecular motor and the Hodgkin and Huxley model [9] for the nervous communications. With the aid of an advanced engineering software (which uses the finite element method to solve di↵erential equations) we will create models that can fully simulate these two processes in order to extract a lot of physical informations about them and to predict a lot of properties that can be, in future, extremely useful during the design stage of both molecular motors and devices which rely their actions on the nervous communications between active fibres.


References [1] Victor E. Saouma. ”Lecture notes in Structural Engineering” (PDF). University of Colorado. Retrieved 2007-11-02. [2] Galileo Galilei, Dialogues Concerning Two New Sciences. Translated from the Italian and Latin into English by Henry Crew and Alfonso de Salvio. With an Introduction by Antonio Favaro (New York: Macmillan, 1914). [3] Chapman, Allan.

England’s Leonardo: Robert Hooke and the Seventeenth-

century Scientific Revolution, (2004), Institute of physics publishing. [4] I. Newton, The mathematical principles of natural philosophy, (1687), Translated into English by A. Motte, 1729. Reprinted in 1968 by Dawsons of Pall Mall, London [5] Engineers’ Council for Professional Development. (1947). Canons of ethics for engineers [6] [7] T. Clark, A Handbook of Computational Chemistry, (1985), Wiley, New York. [8] Risken H., The Fokker-Planck Equation: Methods of Solutions and Applications, 2nd edition, Springer Series in Synergetics, Springer, ISBN 3-540-61530-X [9] Hodgkin, A., Huxley, A. A quantitative description of membrane current and its application to conduction and excitation in nerve. (1952). The Journal of Physiology, 500–544.




Partial Di↵erential Equations: applications and solution Modern science is based on the application of mathematics. It is central to modern society, underpins scientific and industrial research, and is key to our economy. Mathematics is the engine of science and engineering and provides the theoretical framework for biosciences, for statistics and data analysis, as well as for computer science. New discoveries within mathematics a↵ect not only science, but also our general understanding of the world we live in. Problems in biological sciences, in physics, chemistry, engineering, and in computational science are using increasingly sophisticated mathematical techniques. In this chapter we will briefly explain what partial di↵erential equations are, why obtain solutions of them is so complicated and how this problem can be overcomed with appropriate numerical methods.

Partial Di↵erential Equations; who and why Many natural, human or biological, chemical, mechanical, economical or financial systems and processes can be described at a macroscopic level by a set of partial di↵erential equations (i.e. PDEs) governing averaged quantities such as density, temperature, concentration, velocity. Most models based on PDEs used in practice have been introduced in science in the first decades of the nineteenth century in order to study gravitational and electric


fields and to model di↵usion processes in Physics [1]. All the PDEs, in spite of their complexity, involves the first and second partial derivatives. Usually, first derivative involves time, irreversibility of a process, and the second deals with space highlighting symmetries. Nonetheless, PDE theory is not restricted to the analysis of equations of two independent variables and interesting equations are often nonlinear. However, these equations, which were primarily created to model physical processes, played an important role in almost all branches of mathematics and, as a matter of fact, can be viewed as a chapter of applied mathematics [2, 3]. A di↵erential equation is an equation for an unknown function of several independent variables (or functions of these variables) that relates the value of the function and of its derivatives of di↵erent orders. A di↵erential equation is called ordinary di↵erential equation (ODE) if the function u depends only on a single independent variable. A partial di↵erential equation, PDE, is a di↵erential equation in which the unknown function u is a function of multiple independent variables and of their partial derivatives. We can also define as equilibrium equation an equation which models a system evolving with respect only to the space variable. A model can be defined dynamical if the time variable appears; in this case the equation will evolve not only in the space but also in time. Both time and space coordinates are usually independent variables. We can say that Partial di↵erential equations involve rate of change with respect to continuous variables. For example the position of a rigid body is specified by six numbers (i.e. degrees of freedom), but the configuration of a fluid is given by the continuous distribution of several parameters, such as the temperature and pressure. The dynamic for the rigid body problem take place in a finite-dimensional configuration space; the dynamics for the fluid occur in an infinite-dimensional configuration space. This distinction usually makes PDEs much harder to solve than ODEs but, in the case of linear problems, the solution can be simply obtained. The order of a di↵erential equation is that of the highest order derivative that appears in equation; the majority of di↵erential equations arising in applications, in science, in engineering, and within mathematics itself, are of either first or second order, with the latter being by far the most prevalent. Third order equations arise when modelling waves in dispersive media (for example water waves or plasma waves).


Fourth order equations show up in elasticity, particularly plate and beam mechanics. Equations of order

are very rare.

A relatively simple PDE is @u (x, y) = 0 @x


This relation implies that the function u(x, y) is independent of x. However, the equation gives no information on the function’s dependence on the variable y. Hence the general solution of this equation is u(x, y) = f (y)


where f in an arbitrary function of y. The analogous ordinary di↵erential equation (the so called ODE) is du (x) = 0 dx


u(x) = c


which has the solution

where c is any constant value. These two examples illustrates that general solution of ODEs involve arbitrary constants, but solutions of PDEs involves arbitrary functions; a solution of PDE is generally not unique. By a solution of a di↵erential equation we mean a sufficiently smooth function u of the independent variables that satisfies the di↵erential equation at every point of its definition domain. Is not necessary required that the solution is defined for all the possible values of the independent variables but, usually, di↵erential equations are imposed on some domain D contained in the space of independent variables and the needed solution is defined only on D. In general, the domain D will be an open subset usually connected and often bounded with a boundary, denoted @D. A function is called smooth if it can be di↵erentiated sufficiently often, at least so that all of the derivatives appearing in the equation be well-defined on the domain of


interest D. More specifically, if the di↵erential equation has order n, then is required that the solution u be of class C n , which means that it and all its derivatives of order  n are continuous functions in D, and such that the di↵erential equation that relates the derivatives of u holds throughout D.

Boundary conditions and Initial Conditions The general solution to a single nth order ordinary di↵erential equation depends on n arbitrary constants. The solutions to a partial di↵erential equations are yet more numerous, in that they depend on, as said before, arbitrary functions. The solutions to dynamical ordinary di↵erential equations are singled out by the imposition of initial conditions, resulting in an initial value problem. On the other hand, equations which models equilibrium phenomena require boundary conditions to uniquely specify their solutions, resulting in a boundary value problem. A similar specification of auxiliary conditions applies to partial di↵erential equations. Equations modelling equilibrium phenomena are supplemented by boundary conditions imposed on the boundary of the domain of interest. In favourable circumstances, the boundary conditions serve to single out a unique solution. For example the equilibrium temperature of a body is uniquely specified by its boundary behaviour. If the domain is unbounded one must also restrict the nature of the solution at large distances by asking that it remain bounded. Also in the case of PDEs this kind of problem is called boundary value problem. There are three principal types of boundary value problems that arise in most applications [4]. Specifying the value of the solution along the boundary of the domain is called a Dirichlet boundary condition, to honour the nineteenth century analyst Johann Peter Gustav Lejeune Dirichlet. Specifying the normal derivative of the solution along the boundary results in a Neumann boundary condition, named after his contemporary Carl Gottfried Neumann. For example, in thermal equilibrium, the Dirichlet boundary value problem specifies the temperature on its boundary, and our task is to find the interior temperature distribution by solving an appropriate PDE. In the same way, the Neuman boundary value problem prescribes the heat flux through the boundary. In particular, an insulated boundary has no heat flux, and hence the normal derivative of the temperature is zero on the boundary. Prescribing a function along a boundary


and the normal derivative along the remainder will result in a mixed boundary value problem; this kind of boundary conditions in the example of the heat flux will prescribe the temperature along part of the boundary and the heat flux along the remainder. For partial di↵erential equations modelling dynamical process, in which the time is one of the independent variables, the solution is to be specified by one or more initial conditions. The number of initial conditions required depends on the highest order time derivative that appears in the equation. For example, in thermodynamics, which only involves the first order time derivative of the temperature, the initial condition requires to specify the initial temperature of the body at the initial time. Newtonian mechanics describes the acceleration (which is second time derivative of the motion) and so requires two initial conditions; the initial position and the initial velocity of the system. Also in the case of dynamical processes, in bounded systems, one must also impose suitable boundary conditions in order to achieve an unique solution and hence describe the subsequent dynamical behaviour of the physical system.

Numerical methods for solving di↵erential equations Solve analytically di↵erential equations is possible but often difficult [5]. However this is true only for a small number of equations and PDEs that can be solved by explicit analytic formulae are few and far between. Consequently, the development of accurate numerical approximation schemes is essential for extracting quantitative informations as well as achieving a qualitative understanding of the various behaviours of their solutions. In all numerical solutions the continuous partial di↵erential equation (PDE) is replaced with a discrete approximation; discrete means that the numerical solution is known only at a finite number of points in the physical domain. The number of those points can be selected by the user of the numerical method and, in general, increasing the number of points not only increases the resolution (i.e., detail), but also the accuracy of the numerical solution.


The mesh is the set of locations where the discrete solution is computed. These points are called nodes, and if one were to draw lines between adjacent nodes in the domain the resulting image would resemble a net or mesh. Two key parameters of the mesh are

x, the local distance between adjacent points in space, and

distance between adjacent time steps. For the simplicity, lets consider

t, the local x and

t uni-

form throughout the mesh. Many of the basic numerical solution schemes for partial di↵erential equations can be fit into two broad methods. The first is the finite di↵erence method, obtained by replacing the derivatives in the equation by appropriate numerical di↵erentiation formulae. The second category of numerical solution techniques is the finite element method, which will be the method used by us and which will be explained after a brief description of the finite di↵erences method.

Finite Di↵erences method The finite di↵erence method is one of several techniques for obtaining numerical solutions to Equation [6]. The core idea of the finite-di↵erence method is to replace continuous derivatives with di↵erence formulas that involve only the discrete values associated with positions on the mesh. Applying the finite-di↵erence method to a di↵erential equation involves replacing all derivatives with di↵erence formulas. In the heat equation, for example, there are derivatives with respect to time, and derivatives with respect to space. Using di↵erent combinations of mesh points in the di↵erence formulas results in di↵erent schemes. In the limit as the mesh spacing ( x and

t) goes to zero, the numerical solution ob-

tained with any useful scheme will approach the true solution to the original di↵erential equation. In general, a finite di↵erence approximate the value of some derivative of a scalar function u(x0 ) at a point x0 , say u0 (xo ) or u00 (x0 ). on a suitable combination of sampled functions values at nearby points. The underlying formalism used to construct these approximation formulae is known as the calculus of finite di↵erences.


The simplest finite di↵erence approximation is the ordinary di↵erence quotient. u(x + h) h


⇡ u0 (x)


that appears in the original calculus definition of the derivative. Indeed, if u is di↵erentiable at the point x, then u(x) is, by definition, the limit, as h ! 0 of the finite

di↵erence quotients. Geometrically, the di↵erence quotient measures the slope of the secant line through the two points (x, u(x)) and (x + h, x(x + h)) on its graph. For small enough h, this should be a reasonably good approximation to the slope of the tangent line, u0 (x). The step size, h, can be either positive or negative and is assumed to be very small h << 1. When h > 0, is referred to as a forward di↵erence, while u < 0 yields a backward di↵erence.

Finite Element method In the previous sections we have seen the oldest and, in many cases simplest, class of numerical algorithms for approximating the solutions to PDEs, the Finite Di↵erences Method. In the present section we are going to introduce the second of the major numerical algorithms: the FEM, Finite Element Method [7]. The Finite Element Method (FEM) is probably the most powerful method known for the numerical solution of boundary- and initial-value problems characterized by partial di↵erential equations. Consequently, it has a monumental impact on virtually all areas of engineering and applied science. This solution method rely on a more sophisticated understanding of the partial di↵erential equation, in that, unlike finite di↵erences, they are not obtained by simply replacing derivatives by their numerical approximations. Rather, they are initially founded on an associated minimization principle which characterizes the unique solution to a positive definite boundary value problem. The basic idea is to restrict the minimizing functional to an appropriately chosen finite-dimensional subspace of functions. Such a restriction produces a finite-dimensional minimization problem, that can then be solved by numerical linear algebra. When properly formulated, the restricted finite-dimensional minimization problem will have a solution that well approximates the true minimizer, and hence the solution to the original boundary value problem.


There are two fundamental attributes of the method that are at the heart of its great utility and success. Firstly, it is based on the idea of partitioning bounded domains ⌦ into a number N of small, non-overlapping subdomains, the finite elements, over which functions are approximated by local functions, generally polynomials. Secondly, the boundary- and initial-value problems, to which the method is applied, are formulated in an integral form (also called weak form), so that the contributions of each subdomain to the global integrals sum up to produce an integral characterizing the problem over the whole domain. The subdivision of the whole domain into finite simpler sub-domains has several advantages [8]; first of all it permites an accurate representation of complex, multidimensional, geometries; it can simulates anisotropic material properties, it can easy represent the total solution and, finally, with it is possible to capture local e↵ects. In a typical work out of this method one has firstly to divide the domain of the problem into a collection of subdomains, with each of them is represented by a set of element equations to the original problem; the element equations are simple equations that locally approximate the original complex equations to be studied, where the original equations are often partial di↵erential equations (PDEs). This approximation process, in mathematics language, consists in construct an integral of the inner product of the residual and the weight functions and set this integral to zero. In simpler terms, it is a procedure that minimizes the error of approximation by fitting trial functions into the PDE. The residual is the error caused by the trial functions, and the weight functions are polynomial approximation functions that project the residual. The process eliminates all the spatial derivatives from the PDE, thus approximating the PDE locally with a set of algebraic equations for steady state problems and, for transient problems, with a set of ordinary di↵erential equations. These equation sets are the so-called element equations. They are linear if the underlying PDE is linear, and vice versa. Algebraic equation sets that arise in the steady state problems are solved using numerical linear algebra methods, while ordinary differential equation sets that arise in the transient problems are solved by numerical integration using standard techniques such as Euler’s method or the Runge-Kutta method.


Once defined the element equations, the second step to follow in order to solve a PDE with the aid of the FEM is to systematically recombining all sets of element equations into a global system of equations for the final calculation. The global system of equations has known solution techniques, and can be calculated from the initial values of the original problem to obtain a numerical answer. In this step, the global system of equations is generated from the element equations through a transformation of coordinates from the subdomains local nodes to the domain global nodes. This spatial transformation includes appropriate orientation adjustments as applied in relation to the reference coordinate system. The process is often carried out by FEM software using coordinate data generated from the subdomains. FEM can be better understood from its practical application, known as finite element analysis (FEA) [9]. FEA is applied in engineering as a computational tool for performing engineering analysis. It includes the use of mesh generation techniques for dividing a complex problem into small elements, as well as the use of software program coded with FEM algorithm. In applying FEA, the complex problem is usually a physical system with the underlying physics such as the Euler-Bernoulli beam equation, the heat equation, or the Navier-Stokes equations expressed in either PDE or integral equations, while the divided small elements of the complex problem represent di↵erent areas in the physical system. FEA is a good choice for analysing problems over complicated domains (like cars and oil pipelines), when the domain changes (as during a solid state reaction with a moving boundary), when the desired precision varies over the entire domain, or when the solution lacks smoothness. For instance, in a frontal crash simulation it is possible to increase prediction accuracy in ”important” areas like the front of the car and reduce it in its rear reducing in this way the cost of the simulation). Another example would be in numerical weather prediction, where it is more important to have accurate predictions over developing highly nonlinear phenomena (such as tropical cyclones in the atmosphere, or eddies in the ocean) rather than relatively calm areas. Comparison between this method with the finite di↵erence method (FDM) shows the several advantages of the FEM; while FDM in its basic form is restricted to handle rectangular shapes and simple alterations thereof, the handling of geometries in FEM


is theoretically straightforward. The most attractive feature of finite di↵erences is that it can be very easy to implement rather the FEM. Generally, FEM is the method of choice in all types of analysis in structural mechanics (i.e. solving for deformation and stresses in solid bodies or dynamics of structures) while computational fluid dynamics (CFD) tends to use FDM or other methods like finite volume method (FVM). CFD problems usually require discretization of the problem into a large number of cells/gridpoints (millions and more), therefore cost of the solution favours simpler, lower order approximation within each cell. This is especially true for ’external flow’ problems, like air flow around the car or airplane, or weather simulation. Finally, the applications of the FEM; a variety of specializations under the umbrella of the mechanical engineering discipline (such as aeronautical, biomechanical, and automotive industries) commonly use integrated FEM in design and development of their products. Several modern FEM packages include specific components such as thermal, electromagnetic, fluid, and structural working environments. In a structural simulation, FEM helps tremendously in producing sti↵ness and strength visualizations and also in minimizing weight, materials, and costs. FEM allows detailed visualization of where structures bend or twist, and indicates the distribution of stresses and displacements. FEM software provides a wide range of simulation options for controlling the complexity of both modelling and analysis of a system. Similarly, the desired level of accuracy required and associated computational time requirements can be managed simultaneously to address most engineering applications. FEM allows entire designs to be constructed, refined, and optimized before the design is manufactured. This powerful design tool has significantly improved both the standard of engineering designs and the methodology of the design process in many industrial applications [10]. The introduction of FEM has substantially decreased the time to take products from concept to the production line. It is primarily through improved initial prototype designs using FEM that testing and development have been accelerated [11]. In summary, benefits of FEM include increased accuracy, enhanced design and better insight into critical design parameters, virtual prototyping, fewer hardware prototypes, a faster and less expensive design cycle, increased productivity, and increased revenue.


In order to conclude this chapter, we will see which FEA software we have used in order to solve the partial di↵erential equation that underlie the chemical-physics problems we want to solve. It is also of very interest to note that we are able to use, implement and optimize the FEM for describe microscopic problems, very far from the macroscopic, engineering-like, problems for what the FEM have been designed.

COMSOL Multiphysics FEM environment In order to face out our research, as told in the introduction, we want to solve the PDE which describes the dynamic of a molecular motor before and the PDEs which describes the generation and propagation of nervous signals after) we have used the COMSOL Multiphysics [12] to help us setting up the di↵erent problems, to solve the equations behind them and, finally to analyse the solutions of them and extract the wanted quantities. COMSOL Multiphysics is a finite element analysis, solver and Simulation software for various physics and engineering applications, especially coupled phenomena (i.e. multiphysics). COMSOL also o↵ers an extensive interface to MATLAB and its toolboxes for a large variety of programming, preprocessing and postprocessing possibilities. In addition to conventional physics-based user interfaces, the real advantage of COMSOL for our purposes is that it allows the user for entering coupled systems of partial di↵erential equations (PDEs). The PDEs can be entered directly using an apposite coefficient form equation (eq 2.6) in which the user can modify each coefficient matching the equation that the solution is needed or using the so-called weak form. ea

@2u @u + da + r( cru 2 @t @t

↵u + ) + au + ru = f


Finally, COMSOL also o↵ers to the user a CAD interface in order to properly draw/design the wanted domain for the equations without any limitations of dimensionality; finally it o↵ers a proper mesh interface for a correct and suitable meshing of the designed domain. In next chapters we will see how we have used the typical engineering approach (PDE+FEM+COMSOL) to analyse, describe and study microscopic phenomena; in


this way we have start a completely new approach, engineering-like, to molecules (and molecular devices) design.

References [1] Sim´eon-Denis Poisson, A Treatise of Mechanics, (1842), Longman and Company, London.

[2] Evans, Lawrence C., Partial Di↵erential Equations, (1949), American Mathematical Society.

[3] Folland, G. B., Introduction to Partial Di↵erential Equations, (1996), 2nd ed. Princeton, NJ: Princeton University Press.

[4] Powers, David L., Boundary Value Problems and partial di↵erential equations, (2006), Elsevier Academic Press.

[5] Kevorkian, J., Partial Di↵erential Equations: Analytical Solution Techniques, (2000), 2nd ed. New York: Springer-Verlag.

[6] Ames, William F., Numerical Methods for Partial Di↵erential Equations, (1992), Academic Press, Inc., Boston, third edition

[7] Reddy, J.N., An Introduction to the Finite Element Method, (2005), Third ed., McGraw-Hill.

[8] source : http : // initeelementmethod

[9] Ramamurty G., Applied Finite Element Analysis, (2010), I.K. International Publishing House Pvt. Ltd.


[10] Hastings, J. K., Juds, M. A., Brauer, J. R., Accuracy and Economy of Finite Element Magnetic Analysis, (1985), 33rd Annual National Relay Conference. [11] McLaren-Mercedes, Vodafone McLaren-Mercedes: Feature - Stress to impress, (2006), Archived from the original on 2006-10-30. Retrieved 2006-10-03 [12] source :




PDE to study the molecular motors: the Fokker-Planck equation In this chapter we will see what natural and artificial molecular motors are and how an appropriate mathematical model, based on the Fokker-Planck equation, can be extremely useful in order to have a qualitative and quantitative description of how these molecule works.

Natural and artificial molecular motors In recent years the scientific community have reserved a lot of attentions on molecular motors; starting from natural molecular motors such as the ATPase and the Kinesin many kind of artificial molecular motors have been developed [1]. Natural molecular motors are assemblies of proteins within the cellular environment of living organisms that, through complex folding and chemical processes, can perform, as real-world macroscopic engines, mechanical motion for various purposes, such as to transport materials or electrical charges within the cytoplasm of a cell or replicate DNA and other compounds. Also muscle contraction and action is carried out by molecular motor proteins (i.e. the Miosin) as the motion of bacteria through a type of propeller-driven swimming motion [2-3]. Most of these natural molecular motors derive chemical energy for motion from the same basic process that organisms use to produce energy for


life support, by the breakdown and synthesis of the compound adenosine triphosphate (ATP). Though on a basic level molecular motors perform many of the same functions as electro-mechanical motors at the macroscopic human scale, the complex nature of protein folding and chemical reactions that a molecular motor relies on to function, required decades of research in order to fully understand their behaviour. Starting from the knowledge on natural molecular motors, in recent years research in nanotechnology at the atomic and molecular scale has focused on taking biological materials and manufacturing artificial molecular motors which can resemble the motors with which everyday engineering is familiar. A prominent example of this was a motor constructed by a team of scientists at the Boston College of Massachusetts in the US in 1999 that consisted of 78 atoms [4], and took four years of work to construct. The motor had a rotating spindle that would take several hours to make one revolution and was designed to rotate in only one direction. The molecular motor relied on ATP synthesis as its energy source and was used as a research platform to understand the fundamentals of transitioning chemical energy into mechanical motion. Similar research has been completed by scientists using carbon to produce synthetic molecular motors powered by light and heat energy, and recent attempts have developed a method for creating a motor that produces a continuous level of rotational torque [5]. In recent years a lot of di↵erent kinds of molecular motor have been synthesized; there are a lot of di↵erent structures with a lot of di↵erent operative mechanisms; in our work we are going to focus our attention on the rotary molecular motors.

Rotary molecular motors Rotary molecular motors are based on cis-trans photoisomerization processes involving, mainly, -N=N-, -C=N- and -C=C- bonds. These processes are ideal in order to obtain light driven operations of molecular machines because they produces structural changes that can be exploited resulting in large amplitude motions in suitably designed molecular systems like rotary molecular motors. These molecules are systems capable of undergoing unidirectional and repetitive rotations under the action of proper external energy inputs. Construction of molecular rotary motors poses several challenges,


particularly because it is difficult to satisfy the unidirectional rotation requirement. The first motor was based [6] on a symmetric biphenianthrylidene; in such a compound, each one of the two helical subunits linked by a double bond can adopt a righthanded or a left-handed helicity with the result of four di↵erent possible stereoisomers. The cis-trans photoisomerization reactions are reversible and occur upon irradiation at appropriate wavelengths. Counterwise, the inversions of helicities, while maintaining a cis or a trans configuration, are not-reversible thermal processes. In this kind of molecule, a sequence of two energetically uphill light-driven isomerization processes and two energetically downhill thermal helix inversion steps are exploited in order to move this molecular rotor unidirectionality. The directionality of rotation is caused by the energetic preference for the methyl substituents next to the stereogenic centre to adopt axial orientation, which is less sterically demanding and so, energetically favourable. More recently, this first rotary molecular motor have been redesigned in order to improve its performances and, specifically, to lower the barriers to the rate-determining thermal steps. With this purpose a second-generation of this kind of motors was developed [7]; in this new kind of systems the unidirectional rotation can be achieved with a single stereogenic centre and it was observed that an increase of the size of atoms X and Y (fig. 3.1) increases the steric crowding in the fjord region, resulting in slower thermal rotation steps. The passage from a six-membered ring to a five-membered ring in the lower half of the molecule (from structure 4 to structure 5 in figure 3.1) caused [9] an increase in rotation speed by a 108 factor. The thermal helix inversion step of the fast molecular motor 5 has an half-life of 3.15 minutes at room temperature. The structure was also modified [10] with the result that motor 6 could be powered by visible light (436nm) instead of UV light. Finally, in motor 7, light driven unidirectional rotation of the rotor unit when the stator is tethered to the surface of gold nanoparticles was recently demonstrated [8]. Furthermore, the motor was employed as a four-state light-triggered switch to construct a molecular gearbox in which the state of the switch controls the thermal rotation of an appended aromatic moiety [11].


Figure 3.1: - Molecular structures of the second generation of photochemical rotary molecular motors. Image credits [7-8].

This kind of molecular motors can be, in future, used as molecular switch in many di↵erent application such as nanocomputers and nanoswitches; Their functioning principle has been incorporated into a prototype nanocar [12]. Finally, the ability of certain second generation rotary molecular motors to act as an asymmetric catalyst has also been demonstrated [13-14]. It is hence very important, and interesting, to fully understand the functioning principles and the physics of these molecules; understanding completely the thermodynamic of the helix inversion process can be in future useful in order to predict the di↵erent properties of the molecules and assist the scientists during the design of these kind of molecules.

Under the hood: The physics of Molecular Motors dynamic How is possible to describe the dynamic of a rotary molecular motor? To answer this question we have to start from the environment in which the molecular motor lives; molecular motors typically operates in solution; every molecule of the solvent, with its thermal fluctuations, collides to the molecules which are in it and exerts a force on them. In the simplest case this force will result in the well known Brownian Motion


phenomenon. In the simplest case the Brownian motion relates particle which are completely free to move and will results in a completely random motion of them. In the case of a molecular motor, the particle are not completely free to move; the system ”lives” inside a potential energy profile; the result is a biased Brownian Motion and in this case, a molecular motor can also be named as brownian motor. In this phenomena the movement of the solvent molecules will drive the molecular motor in his movement toward equilibrium. In order to correctly model these systems we need a proper physical background which has to start from the definition of Brownian motion.

Brownian Motion Brownian motion is the random motion of particles suspended in a fluid resulting from their collision with the quickly moving atoms or molecules composing it (the fluid). The term ”Brownian motion” can also refer to the mathematical model used to describe such random movements, which is often called a particle theory [15]. This transport phenomenon took this name after the botanist Robert Brown, in 1827, while looking through a microscope at particles found in pollen grains in water, noted that the particles moved through the water but was not able to determine the mechanisms that caused this motion. A concrete description of this interesting phenomena arrived in 1905, when Albert Einstein published a paper [16] that explained in precise detail how the motion that Brown had observed was a result of the pollen being moved by individual water molecules thermal fluctuations around their position. This explanation of Brownian motion served as definitive confirmation that atoms and molecules actually exist, and was further verified experimentally by the Nobel prize Jean Perrin [17]. The direction of the force of atomic bombardment is constantly changing, and at di↵erent times the particle is hit more on one side than another, leading to the seemingly random nature of the motion. The Einstein’s theory of the Brownian motion in composed of two main parts; in the first a formulation of a di↵usion equation for Brownian particles, in which the di↵usion coefficient is related to the mean squared displacement of a Brownian particle, is made. In the second part the di↵usion coefficient is related to measurable physical quantities.


The first part of Einstein’s argument was to determine how far a Brownian particle travels in a given time interval.[16] Classical mechanics was unable to determine this distance because of the enormous number of bombardments a Brownian particle will undergo, roughly of the order of 1021 collisions per second.[18] Thus Einstein was led to consider the collective motion of Brownian particles; he showed that if ⇢(x, t) is the density of Brownian particles at point x at time t, then ⇢ will obey the Fick’s equation for di↵usion: @⇢ @2⇢ =D 2 @t @x


where D is mass di↵usivity. Assuming that all the particles, at time t0 = 0, starts from the origin, the di↵usion equation has the solution

⇢(x, t) = p

⇢0 e 4⇡Dt

x2 4Dt


With this expression is possible to calculate directly the moments; the first moment is seem to vanish, meaning that the Brownian particle is equally likely to move to the left as it is to move to the right. The second moment is, however, non vanishing, and reads x¯2 = 2Dt


This expression means that the mean squared displacement in terms of the time elapsed and the di↵usivity. From this expression Einstein argued that the displacement of a Brownian particle is not proportional to the elapsed time, but rather to its square root. His argument is based on a conceptual switch from the ensemble of Brownian particles to the single Brownian particle, with this approach we can start to speak in terms of probability distribution of the Brownian particles. Obtained the expression 3.2 the determination of the di↵usion coefficient was central to Einstein’s theory; he determined this quantity in terms of molecular qualities as follows.


The Maxwell-Boltzmann distribution for the configuration of Brownian particles which are in a field of force F (x) (i.e. the gravitational field of earth) reads U

⇢ = ⇢0 e


If the force is constant, as would be true for the gravity which acts on particles, the potential energy is

U = Fx


Now the velocity of a particle that is in equilibrium under the action of the applied force and viscosity is given by the Stoke’s law:


F 6⇡⌘a


We can also derive the current density of the particles, i.e. the number of particles crossing unit area in unit time which is


⇢F 6⇡⌘a


For a normal di↵usion process, the particles cannot be created or destroyed, which means that the flux into one space region must be the sum of particles flux flowing out of the surrounding regions. The can be summarized mathematically by the continuity equation: @⇢ + rJ = 0 @t


Combining eqs. 3.1 with 3.8 and integrate over x will give us


@⇢ ⇢F = @x 6⇡⌘a


and, combining eqs. 3.5 with 3.4 we get 1 @⇢ = ⇢ @x


F kB T


So, comparing equations 3.9 and 3.10 we obtain the definition for the di↵usion coefficient D=

kB T 6⇡⌘a


in which T is the temperature, ⌘ is the viscosity of the liquid and a is the size of the particle. This equation means that large particles would di↵use more gradually than smaller molecules, making them easier to measure. Experimental observation confirmed the numerical accuracy of Einstein’s theory. We understand di↵usion in terms of the movements of the individual particles, and can calculate the di↵usion coefficient of a molecule if we know its size or, viceversa, we can calculate the size of the molecule after experimental determination of the di↵usion coefficient. Thus, Einstein connected the macroscopic process of di↵usion with the microscopic concept of thermal motion of individual molecules. The aforementioned discussion holds true only for pure brownian motion; for a molecular motor, we need to think in terms of biased brownian motion. The molecule of the solvent will acts on the moving parts of the molecular motor as they does for a free Brownian particle but, since the motor works inside an anisotropic energetic environment, this would bias the motion of it caused by the solvent hits. The result would essentially be di↵usion of a particle whose net motion is strongly biased in one direction.

Biased Brownian Motion Before we go deeper into the physics of Brownian motor and to study their basic principles we have to introduce a more general principle that runs Brownian motion: the biased Brownian motion. The basic idea is simple: a system which isn’t in thermal equilibrium tends toward equilibrium. If such a system lives in an asymmetric world, then moving towards equilibrium will usually also involve a movement in space. To keep the system moving, we need to perpetually keep it away from thermal equilibrium, which costs energy - this is the energy that drives the motion. This result is reachable in many ways such as fluctuating temperature, chemical reactions, periodically turning


on and o↵ asymmetric potentials or by forcing periodic a Brownian particle.

Figure 3.2: - This figure illustrates how Brownian particles, which was initially located at the point x0 (lower picture), spreads out when the potential is turned o↵. When the potential is turned on again, most particles are captured again in the attraction point x0 , but also exists a probability to find them in the point x0 + L (hatched area). This will result in a net current to the particles to the right.

All these processes can be described by inducing some external potential U(x) in the system. In the case where our periodic potential U(x) has exactly one minimum and maximum per period L as in figure 3.2, it is quite obvious that if the local minimum is closer to its adjacent maximum to the right (fig. 3.2) a positive particle current, x˙ > 0, will arise. So in devices based on Brownian motion, net transport occurs by a combination of di↵usion and deterministic motion induced by proper applied force fields.

In order to describe these phenomena we start introducing the Langevin theory which will lead us to the stochastic di↵erential equations and, finally, the FokkerPlanck equation. Langevin began by simply writing down the equation of motion for the Brownian particle according to Newton’s law under the assumption that the particle experiences two forces:


(i) a fluctuating force that changes direction and magnitude frequently compared to any other time scale of the system and averages to zero over time; (ii) a viscous drag force that always slows the motions induced by the fluctuation term. This equation, named Langevin motion equation, according to Newton’s second law of motion, is m

d2 x(t) = dt2

dx(t) + F (x, t) + f (t)stoch dt

in which F (x, t) represents some optional external force. The friction term

(3.12) ⇣ x˙ is

assumed to be governed by Stoke’s law which states that the frictional force decelerating a spherical particle of radius a is ⇣ x˙ = 6⇡⌘ax˙


where ⌘ is the viscosity of the surrounding medium. For the fluctuating part f (t)stoch the following assumptions are made: (i) f (t)stoch is independent of the space (x), (ii) f (t)stoch varies extremely rapidly compared to the variation of x(t), (iii) the statistical average over an ensemble of particles, f (t)stoch = 0, since f (t)stoch is so irregular. From solutions to the Langevin equations is possible to retrieve a large number of informations about the described system; however this is fully true in case of purely microscopic considerations. If one wants to to study the properties and the dynamic of bigger systems (systems composed of ⇡N particles) limitations of this model and of its classic approach becomes evident. If we want to fully describe the time evolution

of whole a macroscopic environment we have to solve the Langevin’s motion equation for every particle which is part of such a system and it is also necessary to know, for every particle, initial position and velocity. Doing this is not impossible but requires a huge amount of time and calculation power; is better, in order to describe successfully complex problems, to change approach; is thus convenient to take the helpful and the advantages of the statistical mechanics.


A stochastic approach: the Fokker-Planck equation A more efficient approach to the biased Brownian motion problem is given by the Fokker-Planck equation (FPE ) which is just an equation of motion for the distribution function of fluctuating macroscopic variables. The di↵usion equation 3.1 for the distribution function of an assembly of free Brownian particles is a simple example of such a equation. This equation is useful not only for the description of the Brownian motion but it can be used to explain a lot of di↵erent problems such as mathematics or financial problems [19]; is possible to say that the main use of the Fokker-Planck equation is as an approximate description for any Markov process in which the individual jumps are small. In its simplest and more general form the FPE reads: @ @ @2 P (x, t) = (D1 (x, t)P (x, t)) + 2 (D2 (x, t)P (x, t)) @t @x @x


where P (x, t) indicates the probability density distribution of the stochastic variable being studied (it can be the position or velocity of a population of Brownian particles) and the parameter D1 and D2 , which can be dependent by space, time, both or simple constants, are, respectively, the drift and di↵usion coefficients; they will takes di↵erent forms depending of which is the case being solved. It is also of very interest to notice that, in the special case of a zero drift coefficient the FPE reduces simply to the well known di↵usion equation. Mathematically, the FPE is partial di↵erential equation of the second order, of parabolic kind; it is also knowed as forward Kolmogorov relation. Solutions to the FPE can be achieved analytically only in special cases. A formal analogy of the Fokker–Planck equation with the Schroedinger equation allows the use of advanced operator techniques known from quantum mechanics for its solution in terms of eigenvalues and eigenfunctions for a restricted number of cases [19]. For the Brownian particle problem we can get the value of the drift coefficient out of Langevin equation [20] and the value of the di↵usion coefficient directly from Einstein’s considerations [16]. For the particular situation of a Brownian motor, which is the situation of a biased Brownian motion of a particle under the presence of a potential energy profile is possible to rewrite the FPE obtaining a particular form of it knowed as Klein-Kramer equation. This equation is a special case of the Fokker-Planck


equation and it is a motion equation for a distribution function describing both position and velocity of a Brownian particle under the presence of an external force and it reads. @P (x, v, t) = @t in which

@ @ + @x @v


F (x) m


kB T @ 2 P (x, v, t) m @v 2


is the viscous drag constant, m is the mass of the particle, T is the

solvent temperature, kB is the Boltzmann constant and, finally, F (x) =

mf 0 (x) is

the external force result of the potential energy profile. Is thus possible to a↵ord [21] that the Brownian motor lives in a situation in which the inertial e↵ects are negligible, a situation knowed as low Reynold numbers or high viscous coefficients situation. An equation, particular case of the eq. 3.15, which have been written for this special cases, is the Smoluchowski equation: @P (x, t) 1 = @t m

@ @2 F (x) + kB T 2 P (x, t) @x @x


This last equation, which describes the dynamic, in terms of position, of a probability distribution is the one of interest for us in order to achieve the first target of our research work. From solutions of it is hence possible to retrieve a huge quantity of informations about the dynamic and the physics of the system which is studied as we will see in the next chapter when we apply this mathematical approach to the dynamic of the motion of a rotary molecular motor; we will also apply the most recent computational chemistry aids in order to extract the potential energy path where the motor exploit his action.

References [1] Kay, E. R., Leigh, D. and Zerbetto, F. (2007). Synthetic molecular motors and mechanical machines. Angewandte Chemie (International ed. in English) (Vol. 46, pp. 72–191). [2] M. Schliwa, G. Woehlke, Nature 2003, 422, 759. doi:10.1038/ NATURE01601 [3] Molecular Motors (Ed. M. Schliwa) 2003 (Wiley-VCH: Weinheim).


[4] T. Ross Kelly, Harshani De Silva and Richard A. Silva (1999). Unidirectional rotary motion in a molecular system. Nature 401, 150-152. doi:10.1038/43639;

[5] J. Vicario, M. Walko, A. Meetsma and Ben L. Feringa, Fine Tuning of the Rotary Motion by Structural Modification in Light-Driven Unidirectional Molecular Motors, Journal of the American Chemical Society 2006 128 (15), 5127-5135

[6] N. Koumura, R. W. J. Zijlstra, R. A. van Delden, N. Harada, B. L. Feringa, Nature 1999, 401, 152. doi:10.1038/43646

[7] N. Koumura, E. M. Geertsema, A. Meetsma, B. L. Feringa, J. Am. Chem. Soc. 2002, 124, 5037. doi:10.1021/JA012499I

[8] R. A. van Delden, M. K. J. ter Wiel, M. M. Pollard, J. Vicario, N. Koumura, B. L. Feringa, Nature 2005, 437, 1337. doi: 10.1038/NATURE04127

[9] J. Vicario, A. Meetsma, B. L. Feringa, Chem. Commun. 2005, 5910. doi:10.1039/B507264F

[10] R. A. van Delden, N. Koumura, A. Schoevaars, A. Meetsma, B. L. Feringa, Org. Biomol. Chem. 2003, 1, 33. doi:10.1039/B209378B

[11] M. K. J. ter Wiel, R. A. van Delden, A. Meetsma, B. L. Feringa, Org. Biomol. Chem. 2005, 3, 4071. doi:10.1039/B510641A

[12] En Route to a Motorized Nanocar Jean-Fran¸cois Morin, Yasuhiro Shirai, and James M. Tour Org. Lett.; 2006, 8, 1713-1716.

[13] Dynamic Control of Chiral Space in a Catalytic Asymmetric Reaction Using a Molecular Motor Science 18 March 2011: Vol. 331 no. 6023 pp. 1429-1432 doi:10.1126/science.1199844


[14] Heat and Light Switch a Chiral Catalyst and Its Products Science 18 March 2011: Vol. 331 no. 6023 pp. 1395-1396 doi:10.1126/science.1203272 [15] M¨orters, Peter; Peres, Yuval (25 May 2008). Brownian Motion ¨ [16] Einstein, Albert (1905). ”Uber die von der molekularkinetischen Theorie der W¨ arme geforderte Bewegung von in ruhenden Fl¨ ussigkeiten suspendierten Teilchen”. Annalen der Physik 17 (8): 549–560. Bibcode:1905AnP...322..549E. doi:10.1002/andp.19053220806. [17] Perrin J. B. (1916) Atoms [18] Chandresekhar, S. Stochastic problems in physics and astronomy. Rev. Modern Phys. 15, (1943). 1–89. [19] Risken H., The Fokker-Planck Equation: Methods of Solutions and Applications, 2nd edition, Springer Series in Synergetics, Springer, ISBN 3-540-61530-X [20] W.T. Co↵ey, Yu. P. Kalmykov, J. T. Waldron, The Langevin equation, World Scientific, Singapore, 1998 [21] E. M. Purcell, Life at low Reynolds number



Rotary molecular motors dynamic revealed In the previous chapter we have described what the molecular motors are and how their physics can be explained using the Fokker-Planck equation. In this section the first main target of our work is presented; we will provide a general framework that makes possible the estimation of time-dependent properties of a stochastic system moving far from equilibrium like the motion of an artificial molecular motor. The process is investigated and discussed in general terms of non-equilibrium thermodynamics. The approach is simple and can be exploited to gain insight into the dynamics of any molecular-level machine. As a case study, we examine the dynamics of an artificial molecular rotary motor, similar to the inversion of a helix, which drives the rotor from a metastable state to equilibrium as we seen in chapter 3. The energy path that the motor walks was obtained from the results of atomistic calculations. The changes in time of the motor entropy, internal energy, free energy, net-exerted force by the motor, and the efficiency of its unidirectional motion are given, starting from the solution of the Smoulochowski’s equation (solutions obtained with the FEM computational method). The amount of available energy converted to heat due to the net-motion is rather low revealing that the motion is mainly subject to randomness.


Theoretical framework As mentioned in chapter 3 the dynamic of a system (an artificial molecular motor) evolving from a far from equilibrium situation and subject to the environmental fluctuations that acts on it will be done using the Smoluchowski equation @P (x, x0 ; t) @ 2 P (x, x0 ; t) 1 @ =D + @t @x2 @x

@V (x) (P (x, x0 ; t)) @x


where P (x, x0 ; t) is the probability density function related to the position of a particle population, D is the di↵usion coefficient, V (x) is the potential energy profile of the path where the system lives, and

is the friction coefficient of the solvent. According

to Stoke’s law, for spherical objects

= 2d⇡⌘r, with ⌘ the temperature dependent

viscosity coefficient of the medium, r the radius of the rotor, and d the dimensionality of the motion. Combining the continuity equation (4.2) @P (x, x0 ; t) @J(x, t) + =0 @t @x


with equation 4.1, the resulting flux density (or probability current) in units of (s



J(x, t) =


@P (x, x0 ; t) @x

1 @V (x) P (x, x0 ; t) @x


According to Gibb’s definition, the entropy production of the rotor is

S(t) =




P (x, x0 ; t)ln(P (x, x0 ; t))dx



The time derivative of the eq. 4.4 in terms of the probability density, the flux, and the conservative force, F (x) = @x V (x), reads [1] dS = kB dt


L 0

J 2 (x, t) dx DP (x, x0 ; t)

kB D



F (x)J(x, t)dx



The second part at the right hand side of eq. 4.5,

kB D

RL 0

F (x)J(x, t)dx, multiplied

by the temperature of the bath is the heat dissipation rate [2]. It relates to the rate of


change of the entropy of the medium. The entire entropy production rate (entropy of the rotor plus environment) is [1] dStot = kB dt


L 0

J 2 (x, t) dx DP (x, x0 ; t)



The inequality of eq. 4.6 becomes an equality only when the system is at equilibrium and the total flux is zero (detailed balance). In terms of the probability density function the change of the internal energy reads: U=


L 0

V (x) {(P (x, x0 ; t)

P (x0 ; 0)} dx


and the heat dissipated by the particle up to time teq takes the form Qhd =

kB T D


teq 0



F (x)J(x, t)dxdt



and the change in the free energy, combining eqs.4.4 and 4.7, reads

F =


L 0

V (x) {(P (x, x0 ; t) P (x0 ; 0)} dx+ Z L kB T P (x, x0 ; t)ln(P (x, x0 ; t))dx

kB T





P (x0 , 0)ln(P (x0 , 0))dx 0

where teq is the time needed to the system to reach equilibrium. This time can be provided by the analysis of the mean displacement, which connects the initial and the final position of the particle at time t.
x0 > (t) =




x0 )P (x, x0 ; t)dx



The equilibration time, teq , is reached when the mean displacement reacher a plateau and its value is ”almost equal” to the distance between the two minima. The time derivative of the mean displacement provides additional information, namely, the mean velocity of the motion of the motor v(t) =

x0 > (t) = dt



J(x, t)dx



When the mean velocity di↵ers from zero, a net-motion occurs and part of the available energy is converted into heat. The work done by the motor up to time t is given as


the path integral of the net-fore exerted by the motor over the trajectory followed by the R x(t ) system, W = x(0)eq Fnet (x(t))dx(t), where we approximate Fnet (t) = Ff r (t) = v(t). Since the system evolves from a situation which is far from equilibrium, the dissipationfluctuation theorem does not hold for this case a di↵erent approximation is necessary. Taking into account that two competitive minima are present in the bistable potential and that the force is a vector dependent quantity, we write for the net-exerted force by the rotor due to the ith minimum [3] i Fnet (t)


Z t (Z 0



F (x)p(x, t )dx b




F (x)p(x, t )dx dt0




in which p(x, t) is the probability distribution and NOT the probability density function, b is the position of the boundary (to the left or to the right of the minimum), xi,m is the position of the minimum, and xb is the position of the energy barrier. The total net-force exerted by the particle will be deep Fnet (t) = Fnet (t)

shal Fnet (t)


where the superscripts deep and shal indicates the deeper and shallower minimum, respectively. The total work done by the rotor until equilibrium is reached reads W =



Fnet (t)v(t)dt



Coupling eqs. 4.9 and 4.14 we can obtain the efficiency of the motor which, in percent, reads ✏=

W ⇥ 100 F


All the quantities here listed can be easily obtained, for a given potential energy profile and a given solvent, after numerical solutions of eq. 4.1; • the change of the internal energy, eq. 4.7, • the change in the free energy, eq. 4.9, • the mean displacement, eq. 4.10,


• the net exerted force by the particle, eq. 4.12, • the work done by the particle, eq. 4.14, and • the efficiency of the motor, eq. 4.15. In the next sections we will see in detail which is the studied molecular motors, its properties and how we start from experimental evidences in order to construct a proper potential energy profile; we will also see how we solved the Smoluchowski equation which describes the evolution of the system.

The chosen molecular motor The studied molecular motor is described in the article [4]; in that paper a particular kind of rotary molecular motor is synthesized and presented; it is composed by two main parts: a rotor and a stator which are linked by a -C=C- bond which can be photo-isomerized (355nm for 9 seconds) driving the molecule out of equilibrium (unstable situation). The more stable situation will be naturally recovered through a helix inversion process. If we consider that the photo-isomerization occurs instantaneously (

t ⇡ 0), the recovery of the equilibrium will require a time which can varies with the

temperature, the solvent and the substituents presents on the molecule. An example of this kind of molecule is showed in the Figure 4.1.

Figure 4.1: - Scheme of how a rotary molecular motor works. Reprinted with permission from Ref [5]. Copyright 2008 American Chemical Society.

In the latter figure the mechanism that underlie the rotation of the motor is explained; this mechanism is simply to understand talking in terms of potential energy surfaces (PES ); the di↵erent relative positions of the stator and the rotor lies on di↵erent points, which are at di↵erent energies, on a potential energy surface as explained in the Figure 4.2.


Figure 4.2: - Mechanism that underlie the relative rotation of the rotor and the stator in a synthetic rotational molecular motor, Reprinted with permission from Ref [5]. Copyright 2008 American Chemical Society.

For our study we selected the most simple rotary molecular motor which is described in the already cited article; this molecule 4.3 is a cyclopentane based molecular motors in which all the substituents on the six term aromatic carbon rings are hydrogen. Finally there is an S atom at the bottom of the structure.

Figure 4.3: - Structure of the molecule which have been studied in the present work; the red bonds are the ones around which the two parts of the molecule rotates.

The group that has synthesised these molecule have demonstrate that cyclopentanebased molecular rotary motors, which display even less steric hindrance can accomplish


unidirectional rotary motion with a rotational frequencies of the order of the MegaHertz (MHz) at ambient temperatures. Direct measurement of the kinetics of the thermal helix inversion at ambient temperatures is achieved using transient absorption spectroscopy. The experimental data are in full agreement with extrapolated data obtained by analytical techniques and prove that these molecular motors display well-defined rotational behaviour up to MHz frequencies over a wide temperature range. Our computational work starts from the experimental data in order to construct a proper potential energy surface for the molecule; done this, starting from the kinetic data of the molecule the time evolution of the system (evolution on the energy profile) have been reconstructed.

PES construction The kinetic data obtainable experimentally from the molecules can give us a precise idea of the value of the energy barrier which separates the two minima of the potential energy profile and, also of very importance, the velocity in which the system evolves and restores the equilibrium situation. All these data can be obtained directly from the Arrhenius plot showed in Figure and 4.4 which refers to a room temperature transient absorption (i.e. TA) spectroscopy analysis.

Figure 4.4: - Arrhenius plot of the compound in Dichloromethane obtained by TA spectroscopy at ambient temperatures. Reprinted with permission from Ref [4]. Copyright 2008 American Chemical Society.

From these plots is firstly possible to extract the energy barrier activation to thermal helix inversion; for the studied molecule it has a value of 8,4 kcal/mol.


To study of the motion of the molecular motor the determination of a reasonable potential energy profile, which describes as closely as possible the helix inversion process, is required. To reach this goal the Nudged Elastic Band (NEB) method [6] implemented in Amber 11 [7] was used to construct the path of helix inversion process. This computational method requires, as input, the initial and final states (in terms of energetic minima) of the studied transformation (see Figures 4.5 and 4.6), and allow the user to obtain all the intermediate structures and, with the energy of them, to trace the path of minimum energy that joins them; in this way, a total of 100 structures were identified between the initial and final states of the isomerization.

Figure 4.5: - Structure used as a starting point in the calculations.

Figure 4.6: - Structure used as arrive point in the NEB calculations.


The two showed starting points, i.e. the two minima of the potential energy profile, have been obtained modifying the relative position of the stator with respect to the rotor until the energy gains two minima. The first and the last points of the pathway represent a rotation of the rotor (with respect to the stator) of about 43 passing thorough a potential energy maximum.

Once finished the NEB calculations, single point calculations at B3LYP (6-31G(d)) [8] were used to refine the energy of the structures which was obtained. The barrier was 9.5 kcal/mol in good agreement with the experimental value of 8.4 kcal/mole [4]. The bistable potential energy path and its polynomial fit are depicted in the following figure.

Figure 4.7: - Linearised potential energy profile, dots, and polynomial fit, solid line, of the helix inversion path. Three conformational states of the molecule are also shown, one for each minimum and one for the transition state. In red the dihedral angle involved in the rotation.

The angle of the relative rotation, #, varies between 33.58 and

10.55 . The

rotation can be expressed as a function of x = r · # with r = 1nm the radius of

the rotor and the result that x varies from 0 to 4˚ A. With this considerations, the two minima will be located at xA =0.56˚ Aand xB =3.38˚ A, and the total translational motion of the rotor is 2.82˚ A. The potential energy profile was finely fitted by a polynomial of


12th order function that reads:

V (x) =

4 X

ai xj +


3 X

a2i+6 x2j+6



The R-square of this fitting function is 0.999, and the coefficients of the di↵erent orders of the polynomial are listed on table 4.1 ˚n Table 4.1: Polynomial fitting coefficients; the coefficients ai and a2i+6 are given in J/A so that in the polynomial expression distances are in ˚ A and energies in J. a0 5.336 · 10

a1 20

6.995 · 10






1.467 · 1012

6.340 · 1021

a6 6.613 · 1040

a8 6.099 · 1059

The minima of the two well of the system are VA = 1.375 · 10

19 J,



3.107 · 1078

7.3 · 10

20 J

6.471 · 1096

and VB =

while the top of the barrier is set at zero potential energy, VC = 0.

The only parameter that we have to know, in order to go over with the simulations, are the escape rate kAC and the di↵usion and friction coefficients, respectively D and . As told before, the values of the escape rate are taken directly from the experimental evidences and the values of D and

will be parametrized on them starting from real

values of solvent viscosity at the di↵erent temperatures. The Smoluchowski equation, eq. 4.1, is thus solved to reproduce the experimental value of tEV OL , which corresponds to the time needed for 64 % of the probability distribution to be found at the deeper minimum.

Escape rates, friction coefficients and Solution of the Smoluchowski equation As mentioned before, the two coefficients D and

were initially evaluated starting form

real, experimental obtained [9], values of the dichloromethane viscosity. Initial values of these parameters are obtained starting from the fluctuation-dissipation theorem so that D = kB T / where , as told before, is obtained by the well known Stoke’s law. These starting values will be slightly modified to match perfectly the experimental values of the evolution time of the system which can be easily extracted from the Arrhenius plot in fig 4.4. As a consequence, the final di↵usion and friction coefficients


do not obey the dissipation-fluctuation theorem, which is reasonable since the process is not at equilibrium. The results are listed in table 4.2: Table 4.2: Di↵usion and friction coefficients for the molecule in Dichloromethane. Evolution and equilibration times and escape rates are also listed

T (K) 297

D(m2 /s) ⇥ 10 2.495


(N s/m) ⇥ 10


tEV OL (s)

teq (s)

kAC (s






























Once the potential energy profile and the coefficients are ready is possible to continue the study solving the equation related to the evolution of the system. As explained in the chapter 2 the solutions of the FPE have been reached with the COMSOL Multiphysics software which uses the FEM method for solving the di↵erential equations. Once created a proper domain for the equation on which the equation (and the related parameters) are defined, defined proper boundary conditions, and times for the timedependent solutions, the equation is ”simply” solved. Changing the parameters D and ↵ is possible to obtain all the solution for the di↵erent temperature conditions. The very first solution which can be obtained with this method is the space distribution of the probability density (normalized) and every plot represents a successive time in the evolution of the system. At time t0 = 0 the initial probability density function is a Gaussian distribution, positioned at the shallower minimum, that reads

P (x, x0 ; 0) = p where xA = 0, 56˚ A, and

= 97 ⇥ 10

1 e 2⇡

(x xA )2 2 2


13 m.

After a proper numerical analysis of the data given by the FEM environment every wanted property is easy obtained.


Results and analysis of them The very first data which can be extracted is the mean displacement, eq. 4.10, which is showed in figure 4.8a for five di↵erent temperatures. The same information can be also provided by the time evolution of the probabilities at the two minima xA and xB as showed in figure 4.8b.

Figure 4.8: - a) Time evolution of the rotor mean displacement for five di↵erent temperatures. b) Time evolution of the probability at the minimum of each well, respectively A and B for five di↵erent temperatures.

Both the properties give the same trends at all temperatures and indicates equilibration times of the order of 1 µs, the values of teq are listed in table 4.2, and are in agreement with experimental evidences, tEV OL . The mean displacement (MD) reaches a plateau at 2,76 ˚ A, very close to the distance between the two minima, ˚ xB = xA = 2, 82A. Is interesting to note that the value of the plateau is very weakly temperature dependent, while teq is temperature dependent. The higher is the temperature, as expected, the faster is the net-motion of the motor. This finding is also supported by the time evolution of the probabilities at the positions xB and xA where for higher temperatures the equilibration time is shorter. Higher temperatures means greater fluctuations, which helps the motor to pass the energy barrier in shorter times.


The time evolution of the entropy of the rotor is depicted in figure 4.9a, while eq. 4.5 gives the rate of the entropy production of the rotor and is depicted in figure 4.9b.

Figure 4.9: - a) Time evolution of the rotor entropy for five di↵erent temperatures. b) Rate of change of the entropy production of the rotor as a function of time for five temperatures.

At t = 0, the entropy of the rotor is not zero (the initial distribution is very shape but not a dirac function so that the rotor is not fully localized at a specific point with probability 1). For t > 0, the distribution becomes broader exploiting thermal fluctuations and the entropy of the rotor increases, Fig. 4.9a. However, the rate of change of entropy production decreases indicating the role of the energy barrier, which prevents the continuously broadening of the initial distribution. The rate of change has a turning point where the entropy of the rotor takes its maximum value. Such a point corresponds to the time when the two minima have the same probability of occupation, see Fig. 4.9b. When more than 50 % of the initial distribution has moved to the deeper minimum, the entropy of the rotor decreases and reaches a plateau at equilibrium. Thermal fluctuations make the process faster or slower depending on the temperature of the bath, and are furthermore responsible for transduction of energy from the environment to the rotor. This is highlighted on figure 4.9b where the rate of change of entropy is depicted. After the unidirectional motion starts, there is a flow of energy from the environment to the rotor and the rate of change is positive. After the turning point of entropy production, the rotor starts to dissipate energy back to the environment and its entropy decreases. Notice that the mechanism of energy ex-


change between rotor and environment (fluctuations and dissipation) co-exists during the spontaneous process and they cancel each other at equilibrium. The amount of the energy exchanged between the two parts (rotor and environment) is temperature dependent, higher temperatures lead to a grater change of the rate of the entropy of the rotor, which means that there is a faster conversion of the available energy into heat. The process described by Figs 4.8 and 4.9 is irreversible. The strongest criterion to indicate irreversibility is the non-vanishing of the rate of the total entropy production [10]. However, for systems with constant volume and isothermal conditions also the change of the Helmholtz free energy can characterize the process, irreversibility is accompanied by a decrease of the free energy until it reaches a plateau at equilibrium. In Fig. 4.10, the change of the Helmholtz free energy and the change of the internal energy are shown.

Figure 4.10: - a) Time evolution of the Helmholtz free energy for five di↵erent temperatures. b) Mean internal energy as a function of time for five temperatures.

The free energy decreases over time and it remains constant as the system reaches equilibrium indicating the irreversibility of the process, Fig 4.10a. The value of the plateau represents the free energy of the rotor at equilibrium. This value is in agreement with the value obtained from the equation Feq = kB T lnZ + V , where Z = R L V (x)/k T B dx is the partition function, and V = V (xB ) V (xA ). Decreasing the 0 e temperature the changes becomes slower, due to lower thermal fluctuations. the change in internal energy follows the same trends of the change of free energy, Fig. 4.10b.


The net-force exerted by the rotor, eq. 4.12, is depicted in Fig. 4.11.

Figure 4.11: temperatures.

- Net-force exerted by the rotor as a function of time at five di↵erent

The maximum value of the net-force varies between 0.2 and 0.6 pN and depends on the temperature. The higher the temperature the higher the maximum value of the force is. The maximum value of the net-force is approximately ten times smaller than what is usually the case for the nanomachines used by Nature (usually a force of order of some pN).[11] The work done by the motor and its efficiency are listed in Table 4.3. The efficiency of the process is very small; it indicates that the unidirectional rotation is not governed directly by the presence f the two minima, where the second acts as an attractor. Rather it is the result of the redistribution of the initial packet of probability towards the deeper minimum [12-13]. The low efficiency and the low yeald of work produced by the rotor is due to the high entropic losses. Recently, the same trend has been pointed out for a random walker moving between two minima [14].


Table 4.3: Efficiency, maximum work and maximum change of the free energy of the rotor helix inversion for five di↵erent temperatures.

F (J) ⇥ 10


W (J) ⇥ 10

T (K)

✏c (%)






















Conclusion In order to conclude this first part of the thesis is possible to say that, combining wellknown concept of non-equilibrium thermodynamics, which allows the calculation of time dependent properties, with atomistic calculations, which allow the determination of the reaction path, and leveraging on the solution of Smoluchowski’s equation, the probability density function is recorded at all times for an artificial molecular motor. The approach is simple and general and can be applied in the development of artificial molecular motors revealing the role of friction coefficient and the position of the barrier towards reducing randomness.

References [1] H. Qian, Mesoscopic nonequilibrium thermodynamics of single macromolecules and dynamic entropy-energy compensation, Phys. Rev. E 65, 016102, (2001). [2] M-P. Qian, M. Qian, G. L. Gong, The reversibility and the entropy production of Markov processes, Contemp. Math. 118, 255, (1991). [3] E. Bakalis, Designing Nanomachines: A Theoretical and Computational Approach, J. Comput. Theor. Nanosci. 7, 1, (2010). [4] Klok, M., Boyle, N., Pryce, M. T., Meetsma, A., Browne, W. R., & Feringa, B. L.(2008). MHz Unidirectional Rotation of Molecular Rotary Motors, 2, 10484–10485.


[5] Conyard, J, Addison, K, Heisler, IA, Cnossen, A, Browne, WR, Feringa, BL and Meech, SR (2012) Ultrafast dynamics in the power stroke of a molecular rotary motor. Nature Chemistry, 4 (7). pp. 547-551. ISSN 1755-4330 [6] Bergonzo, C., Campbell, A.J., Walker, R.C., Simmerling, C., A Partial Nudged Elastic Band Implementation for Use with Large or Explicitly Solvated Large Systems, Int. J. Quant. Chem. 109, 3781, (2009). [7] David A. Case, T. A. Darden, T. E. Cheatham, Carlos L. Simmerling, J. Wang, Robert E. Duke, Ray Luo, Michael Crowley, Ross C. Walker, W. Zhang, K. M. Merz, B. Wang, S. Hayik, Adrian Roitberg, Gustavo Seabra, I. Kolossvary, K. F. Wong, F. Paesani, J. Vanicek, X. Wu, Scott R. Brozell, Tom Steinbrecher, Holger Gohlke, L. Yang, C. Tan, J. Mongan, V. Hornak, G. Cui, D. H. Mathews, M. G. Seetin, C. Sagui, V. Babin, Peter A. Kollman, University of California, San Francisco [8] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2009 [9] CRC Handbook of Chemistry and Physics, 47th ed. pp F-33 - F-38 [10] I. Prigogine, Introduction to thermodynamics of irreversible processes,, (John Wiley & Sons, 3nd Ed., 1968).


[11] M. E. Fisher, and A. B. Kolomeisky, The Force exerted by a molecular motor, Proc. Natl. Acad. Sci. USA 96, 6597-6602, (1999). [12] M. Kosmas and E. Bakalis, Di↵usive motion in the presence of an array of interacting centers, Phys. Lett. A 358, 354, (2006). [13] E. Bakalis Physica A, Explicit propagators for a random walker and unidirectionality on linear chains, Physica A 391, 3093, (2012). [14] E. Bakalis, F. Zerbetto, Calculating the force exerted by simple molecular machines, MATCH Commun. Math. Comput. Chem. 66, 549, (2011).



PDE and nervous signals: the Hodgkin and Huxley model In the previous sections we have seen how is possible to solve a PDE which describes a physical phenomenon in order to extract important informations about it; in the following sections we want to expand this approach to describe one of the most important and most studied phenomenon in nature: the generation and transmission of the nervous signals between excitable cells. After the description of the nervous cells from the physiological point of view, we will see how is possible to describe their ability to interact with each of them, especially for informations exchange, with a mathematical model, necessarily simplified, based on an electrical analogy between the cells and their properties, representing them in terms of voltages and currents.

Physiology of an excitable cell: an overview Every living organism is composed of cells [1], each surrounded by a membrane that separates it from the surrounding space (i.e. extracellular space). In a living organism there are a lot of di↵erent kinds of cells each of them with di↵erent functions but their membranes have a general common structure: each of them is a very thin film (its thickness varies from 7 to 9 nm); the basic structure of the cell membrane consists


of two phospholipid layers; phospholipids are molecules that in aqueous solution tend to rise to the surface disposing one next to other with the hydrophilic polar ”heads” facing the water and the hydrophobic, non-polar, ”tails” (which are long fatty acids chains) disposed in the opposite direction (see Figure 5.1). The membrane itself also contains a lot of proteins which performs a lot of di↵erent functions such the transport of molecules and ions through the membrane itself. The biologic membranes can be viewed as barriers with the result that the same molecule can have an intracellular concentration di↵erent to the extracellular one.

Figure 5.1: - The cell membrane.

The ion channels can interact through their common transmembrane potentials and capacitances. A remarkable result of these interactions is that such electrical active tissue can generate a transient pulse of electrical changes, i.e. an action potential, across the cell membrane. The action potential cycle consists of a rapid membrane depolarization (i.e., an increase in transmembrane potential) followed by a slower recover to resting conditions. Once an action potential is initiated at one site on an extensive membrane, it initiates action potentials at adjacent sites, thus leading to a sequence of action potentials throughout the remaining membrane. A simple cellular electrophysiological model is that shown in Figure 5.2 [2]. The cell membrane separates the extracellular and intracellular spaces. Both regions may be idealized as passive and uniformly conducting (with di↵erent conductivities). If an appropriate stimulating current is passed between a pair of electrodes across the


membrane of the cell, a remarkable series of events ensues, which may be observed by recording the transmembrane voltage across the membrane as a function of time. The action potential is a consequence of the stimulus, but it is generated by the charged energy stored in the concentration di↵erences that exist across the excitable membrane. The action potential is generated by the membrane’s utilizing this stored energy to allow first the flow of sodium ions (to move the voltage up) and then the flow of potassium ions (to move the voltage down). When one examines this phenomenon in detail, it is seen to consisted of a series of remarkably complex events as we will see in the next sections.

Figure 5.2: - Electrical Stimulation of an excitable cell. The stimulus elicits an action potential through, top of the picture, a stimulator S ; the result of the stimulation is readed by the voltmeter V. The stimulator injects current into the intracellular volume; the voltmeter measures transmembrane voltage as a function of time. Image Credits [2]

Another interesting property of the action potential is its non-linear behaviour; if the stimulus is reduced by half, then no action potential will occurs. Counterwise, if the stimulus is doubled, the action potential deflection remains largely unchanged, but the latency is markedly reduced.


This chapter will consider action potentials from several perspectives; firstly we will analyse the physiological bases of the action potentials generation and propagation [1-11]; then we will move on how the experimental findings leds to the Hodgkin and Huxley mathematical model [17].

Electrical activity of the cell membrane The membrane double layer, which is a very thin electrical insulant layer that separate the extracellular space from the intracellular has, on opposite sides, di↵erent ions with di↵erent charges; the potential di↵erence between inside and outside of the membrane, called membrane potential, is thus determined by the charges that adhere directly to the membrane charging it [1, 3]. A single isolated cell present, at its resting state, a potential di↵erence across the membrane of approximately -70/-100 mV (negative inside); this value of the membrane potential is explainable analyzing the selective permeability of the membrane with respect to the di↵erent ions present in the intra- and extracellular space. The membrane is slightly permeable to the sodium ions Na+ and almost total permeable to the potassium K+ and chloride ions Cl . The total concentration of the potassium ions is much greater inside the cell than outside while, for both the sodium and chloride ions the behaviour is the opposite; their concentration is greater in the extracellular medium than in the intracellular space; the concentration di↵erence is approximately ten-fold for both potassium and sodium ions, so there is about 10 times the potassium ion concentration inside the neuron than outside and there is about 10 times the sodium ion concentration outside the neuron as inside. The big di↵erences in concentration of these two ions results in very strong concentration gradients for both of them which generates a flow (i.e. di↵usion) of ions from higher concentrated region to the less concentrated is generated. This di↵usion phenomena is governed by the Fick’s di↵usion law and depends on the permeability of the membrane (with respect to the single ions) and on the concentration gradient. Such a movement of one ion across the membrane would result in a net imbalance of


charge across the membrane (a membrane potential) and as more ions pass the membrane as the membrane voltage increases. The charge accumulation on the opposite sides of the membrane generates an electric field which induces a charge flow which is opposed to the flow generated by the gradient thus, the building membrane voltage is an increasing force that acts counter to the tendency for net movement of ions down their concentration gradient. The opposite direction of this electric gradient is necessary to maintain the correct concentration of ions on the two sides of the barrier; when the membrane voltage has grown to the extent that its ”strength” matches the concentration gradients the forces which are applied to ions are now the same strength and oriented in opposite directions, the system reaches the equilibrium; the tendency of ions to follow their concentration gradient is matched by the tendency of the membrane voltage to pull them in the opposite direction; the ions will continue to move across the membrane, but the rate at which they enters and leaves the cell are the same, thus, there is no net ionic current. Because the ions are at equilibrium, membrane potential is stable, or ”resting”; in a resting cell the concentration gradient and the electric one are equilibrated and the equilibrium, or resting potential is defined as the value of Vm (the so-called membrane potential) at which there is a zero net current of any ion. The resting voltage, which is the result of the passage of the ions through the membrane, is the result of several ion translocating enzymes located in the plasma membrane, steadily operating in parallel [3]; The two most important types of these enzymes, which are membrane ion transport proteins, are ion channels and ion transporters. Ion channel creates paths across cell membranes through which ions can passively di↵use without direct expenditure of metabolic energy. They have selectivity for certain ions, thus, there are potassium-, chloride-, and sodium-selective ion channels. Ion transporters, also called ion pumps, are transmembrane proteins that moves ions across the membrane against their concentration gradient, in contrast to ion channels, where ions go through passive transport. These primary transporters are enzymes that convert energy from various sources, including ATP, sunlight, and other redox reactions, to potential energy stored in an electrochemical gradient. It is also important to observe that the resting membrane potential is not an equilibrium potential but rather a dynamic potential as it relies on the constant expenditure of energy (for ionic pumps


as mentioned above) for its maintenance.

Membrane electrical properties: Calculating the membrane potential The resting membrane potential is dominated by the ionic species in the system that has the greatest conductance across the membrane. For most cells this is potassium; due to the active transport of potassium ions, the concentration of potassium is higher inside cells than outside. Most cells have potassium-selective ion channel proteins that remain open all the time. There will be net movement of positively charged potassium ions through these potassium channels with a resulting accumulation of excess negative charge inside the cell. The outward movement of positively charged potassium ions is due to random molecular motion (di↵usion) and continues until enough excess negative charge accumulates inside the cell to form a membrane potential which can balance the di↵erence in concentration of potassium between inside and outside the cell. ”Balance” means that the electrical force (potential) that results from the buildup of ionic charge, and which impedes outward di↵usion, increases until it is equal in magnitude but opposite in direction to the tendency for outward di↵usive movement of potassium. This balance point is an equilibrium potential as the net transmembrane flux (or current) of K+ is zero. The equilibrium potential for a given ion depends only upon the concentrations on either side of the membrane and the temperature. It can be calculated using the Nernst equation [4, 5]:

Eeq,K + =

RT [K + ]o ln zF [K + ]i


where Eeq,K + is the equilibrium potential for potassium, measured in Volts; R is the universal gas constant, 8.314 J · K


· mol


T is the absolute temperature expressed in Kelvin;

z is the number of elementary charges of the ion of which the equilibrium potential is wanted F is the Faraday constant, 96.485 C· mol


[K + ]o is the extracellular (subscript o is for outside) of Potassium;


[K + ]i is the intracellular (here i means inside) of Potassium; Potassium equilibrium potentials of around -80 millivolts (inside negative) are common. Di↵erent potentials are observed in di↵erent species, di↵erent tissues within the same animal, and the same tissues under di↵erent environmental conditions. Applying the Nernst Equation above, one may account for these di↵erences by changes in relative K+ concentration or di↵erences in temperature. For common usage the Nernst equation is often given in a simplified form by assuming typical human body temperature (37 C), reducing the constants and switching to Log base 10. For Potassium at normal body temperature one may calculate the equilibrium potential in millivolts as: Eeq,K + = 61.54mV log

[K + ]o [K + ]i


Likewise the equilibrium potential for sodium (Na+ ) at normal human body temperature is calculated using the same simplified constant. The resting potential of a cell can be most thoroughly understood by thinking of it in terms of equilibrium potentials and can be calculated with the Goldman-Hodgkin-Katz voltage equation [6] using the concentrations of ions as for the equilibrium potential while also including the relative permeabilities, or conductances, of each ionic species. Under normal conditions, it is safe to assume that only potassium, sodium N a+ and chloride Cl

ions play large roles for the resting potential: Em

RT = ln F

PN a+ [N a+ ]o + PK + [K + ]o + PCl [Cl ]i PN a+ [N a+ ]i + PK + [K + ]i + PCl [Cl ]o


This equation resembles the Nernst equation but has a term for each ion; note that the extracellular and intracellular concentrations of the chloride ion are inverted as chloride has negative charge. Em is the membrane potential, measured in Volts, R, T and F are, respectively, the universal gas constant, the absolute temperature and the Faraday constant. PX is the relative permeability of ion X and [X]y is the concentration of ion X in the extracellular space (y = o) and in the intracellular space (y = i). Typical value


for the membrane potential of an excitable cell is about -90 mV (for the nervous cell this is approximately -70 mV) and this value can be compared with the equilibrium potentials for the Na, K and Cl ions. Using typical values for the ionic equilibrium potentials (VN a =65mV, VCl =-90mV and VK =-102mV) is possible to note that the total resting potential doesn’t coincide with them; in particular the potassium potential is 12mv more negative than it and this consists in a continue outward potassium flux and, in the same way, a continue inward sodium flow. As mentioned before, all these passive flows are counterbalanced by the ion pumps in the way to maintain constant the intra- and extracellular ionic concentrations.

Membrane electrical properties: The action potential The action potential is a short-lasting event in which the membrane potential of a cell (excitable cells, which include neurons, muscle cells, and endocrine cells, as well as in some plant cells) rapidly rises and falls, following a consistent trajectory (Figure 5.3). In neurons, they play a central role in cell-to-cell communication. In other types of cells, their main function is to activate intracellular processes. In muscle cells, for example, an action potential is the first step in the chain of events leading to contraction. In beta cells of the pancreas, they provoke release of insulin [7].

Figure 5.3: - Action potential phases. Image credits [8]

Action potentials are generated by special types of voltage-gated ion channels embedded in a cell’s membrane. These channels are shut when the membrane potential is


near the resting potential of the cell, but they rapidly begin to open if the membrane potential increases to a precisely defined threshold value. With in mind the Figure 5.3, the di↵erent phases of the generation of the action potentials are [9]: 1 - A stimulus is received by the dendrites of a nerve cell. This causes the Na+ channels to open; this allow an inward flow of sodium ions, which changes the electrochemical gradient and if the opening is sufficient to drive the interior potential from -70 mV up to -55 mV (the so-called threshold potential), the process continues. 2 - Having reached the threshold, more Na+ channels (also called voltage-gated channels) open. The Na+ influx drives the interior of the cell membrane up to about +30 mV. The process to this point is called depolarization. The process then proceeds explosively until all of the available ion channels are open, resulting in a large upswing in the membrane potential. 3 - The rapid influx of sodium ions causes the polarity of the membrane to reverse, and the ion channels then rapidly inactivate. As the sodium channels close, sodium ions can no longer enter the neuron, and they are actively transported out of the plasma membrane. At this point the K+ channels open. Since the K+ channels are much slower to open, the depolarization has time to be completed. Having both Na+ and K+ channels open at the same time would drive the system toward neutrality and prevent the creation of the action potential. 4 - With the K+ channels open, the membrane begins to repolarize back toward its rest potential. At this point there is an outward current of potassium ions, returning the electrochemical gradient to the resting state. 5 - The repolarization typically overshoots the rest potential to about -90 mV. This is called hyperpolarization and would seem to be counterproductive, but it is actually important in the transmission of information. Hyperpolarization prevents the neuron from receiving another stimulus during this time, or at least raises the threshold for any new stimulus. Part of the importance of hyperpolarization is in preventing any stimulus already sent up an axon from triggering another action potential in the opposite direction. In other words, hyperpolarization assures that the signal is proceeding in one direction.


6 - After hyperpolarization, the Na+ /K+ + pump eventually brings the membrane back to its resting state of -70 mV. In animal cells, there are two primary types of action potentials, one type generated by voltage-gated sodium channels, the other by voltage-gated calcium channels. Sodium-based action potentials usually last for under one millisecond, whereas calciumbased action potentials may last for 100ms or longer. In some types of neurons, slow calcium spikes provide the driving force for a long burst of rapidly emitted sodium spikes. In cardiac muscle cells, on the other hand, an initial fast sodium spike provides a ”primer” to provoke the rapid onset of a calcium spike, which then produces muscle contraction [1, 3].

Membrane electrical properties: Action potential propagation Once an action potential is generated (by an opportune stimulus) at one membrane site, the action potential itself will be the stimulus for a neighbouring segment of the cell membrane. When the propagating (over the membrane surface) signal reaches the axon, it proceeds down this ”transmission line” by successive excitation of segments of the axon membrane [2].

Figure 5.4: - Scheme of the propagation of a nervous signal along an axon. Image credits [8].

Just the successive stimulation of action potentials would result in slow signal transmission down the axon. In order to enable fast and efficient transduction of electrical signals in the nervous system, certain neuronal axons are covered with myelin sheaths [1]. Myelin is a multilamellar membrane that enwraps the axon in segments separated


by intervals known as nodes of Ranvier. The propagation speed is considerably increased by the action of these sheath which prevents the gates on that part of the axon from opening and exchanging their ions with the outside environment. Myelin sheath reduces membrane capacitance and increases membrane resistance in the inter-node intervals, thus allowing a fast movement of the signal. At the uncovered areas of the axon membrane (i.e. the nodes of Ranvier), the ion exchange necessary for the production of an action potential can take place. The action potential at one node is sufficient to excite a response at the next node, so the nerve signal can propagate faster by these discrete jumps than by the continuous propagation of depolarization/repolarization along the membrane. This enhanced signal transmission is called saltatory conduction. [18]

Figure 5.5: - Action potential propagation along a myelinated axon. Image credits [8].

The axon is made up of connected segments of length about 2mm and its diameter is typically 20 µm. Axon diameters may vary from 0.1 µm to 20 µm and may be up to a meter long. The much-studied squid has a giant axon of about a millimeter in diameter [10]. The action potential travels along the axon at speeds from 1 to 100 m/s and, in general, increases with axonal diameter; for axons larger than a minimum diameter (roughly 1 µm), myelination increases the conduction velocity of an action potential, typically ten times. Conversely, for a given conduction velocity, myelinated fibers are smaller than their unmyelinated counterparts. For example, action potentials move approximately at the same speed (25m/s) in a myelinated frog axon and in an unmyelinated squid giant axon, with but the frog axon 30-fold smaller in diameter and 1000-fold smaller in cross-sectional area. Also, since the ionic currents are confined to the nodes of Ranvier, far fewer ions ”leak” across the membrane, saving metabolic


energy. The myelin sheaths are about 1mm in length and also the length of myelinated segments is important to the success of saltatory conduction. They should be as long as possible to maximize the speed of conduction, but not so long that the arriving signal is too weak to provoke an action potential at the next node of Ranvier. In nature, myelinated segments are generally long enough for the passively propagated signal to travel for at least two nodes while retaining enough amplitude to fire an action potential at the second or third node. Some diseases degrade myelin and impair saltatory conduction, reducing the conduction velocity of action potentials. The most well-known of these is multiple sclerosis, in which the breakdown of myelin impairs coordinated movement.

Membrane electrical properties: Experimental methods for measuring action potentials The study of action potentials has required the development of new experimental methods. The initial work, prior to 1955, focused on three goals: isolating signals from single neurons or axons, developing fast, sensitive electronics, and shrinking electrodes enough that the voltage inside a single cell could be recorded. The first problem was solved by studying the giant axons found in the neurons of the squid genus Loligo [11, 12]; these axons are so large in diameter (approx. 1mm, 100-fold larger than a typical neuron) that they can be seen with the naked eye, making them easy to extract and manipulate. However, the Loligo axons are not representative of all excitable cells, and numerous other systems with action potentials have been studied. The second problem was addressed with the crucial development of the voltage clamp,[11] which allows scientists to study the ionic currents underlying an action potential in isolation, and eliminated a key source of electronic noise, the capacitive current associated with the capacitance Cm of the membrane. Since the current equals Cm times the rate of change of the transmembrane voltage Vm, the solution was to design a circuit that kept Vm fixed (dV /dt = 0) regardless of the currents flowing across the membrane. Thus, the current required to keep Vm at a fixed value is a direct reflection of the current flowing through the membrane.


The third problem, that of obtaining electrodes small enough to record voltages within a single axon without perturbing it, was solved in 1949 with the invention of the glass micropipette electrode; with this method electrode tips that are as fine as 10nm, which also confers high input impedance, can be produced. While glass micropipette electrodes measure the sum of the currents passing through many ion channels, studying the electrical properties of a single ion channel became possible only in the 1970s with the development of the patch clamp by Erwin Neher and Bert Sakmann [13]. For this they were awarded the Nobel Prize in Physiology or Medicine in 1991 [14]. Patchclamping verified that ionic channels have discrete states of conductance, such as open, closed and inactivated.

Voltage Clamp The total membrane potential, at resting and during an action potential, is the result of ionic movements through the membrane. These movements generates a current, i.e. membrane current which is the sum of each individual ionic current. In order to describe the activation mechanism quantitatively, one must be able to measure selectively the flow of each constituent on of the total membrane current especially the sodium and potassium current. The voltage clamp (and also space clamp but we will discuss it later) were techniques introduced by Hodgkin and Huxley and were crucial to the task of separating the capacitive, sodium and potassium components of membrane current in their measurements on the giant squid axon [15]. The voltage clamp is a feedback arrangement where the transmembrane potential is held constant electronically during an action potential or subthreshold response. It was conceived as a way to eliminate the complication of the capacitive current since, if dVm /dt = 0, then, obviously, Cm dVm /dt = 0 and the capacitive current is zero. However, even with the voltage clamp separation of ionic fluxes, the quantitative description of the total membrane voltage through each ionic constituent required (by Hodgkin and Huxley) a particular application of the Nernst-Planck equation which led to the famous HH model which can describe quantitatively both the active and passive


membrane behaviour and to understand the membrane mechanism that underlie its electrophysiological behaviour.

Target: separation into individual ionic components The voltage clamp was carefully designed with particular goals in mind, and made possible by a carefully constructed experimental apparatus [12, 16] which was highly successful in providing the needed data. These data allowed the investigators to separate the action potential, a composite event, into component currents, which provided the basis for understanding what was actually happening during the action potentials time course. In the earlier era of the neurophysiology, figuring out how to determine the individual components was the challenge faced by Hodgkin, Huxley, and other investigators of their time. The components of the transmembrane current during an action potential includes the ionic flux plus a capacitive current. Since the capacitance is fixed during an action potential, its current IC , is given by IC =d(Cm Vm )/dt=Cm dVm /dt


Consequently, in a circuit arranged to apply a voltage step across the entire membrane, the capacitive current component will be absent after the voltage step. The removal of capacitive current simplifies the analysis of the remaining current components, because they must then consist entirely of ionic components. Hodgkin and Huxley (HH) reasoned successfully, based on experiments described earlier, that the chloride contribution to the total current did not need to be included explicitly. The major task that remained was to separate the ionic flux into its sodium and potassium components. This separation turns out also to be facilitated by the measurements of current under constant transmembrane potential conditions. To this end, the


voltage-clamp and space-clamp device illustrated in the next figure was developed. In the voltage clamp as designed by Hodgkin and Huxley, a simple proportional controller (as shown in Figure 5.6) is used to keep the membrane potential at a preset value.

Figure 5.6: - Schematic diagram showing the voltage and space-clamp apparatus as developed by HH. Image Credits [16].

Voltage Clamp measurements A typical record resulting from the application of a step change in membrane voltage is shown in Figure 5.7 [2]. In the lower panel of the figure, one notes an early inward current followed by a rise to an asymptotic outward current. An initial capacitive surge is completed in 20 µsec, corresponding to the presence of a capacitor with C = 1.0 µF/cm. Because of the very short time constant, this current drops to zero before the ionic current becomes significant, and hence is normally ignored in studies of the latter. The initial flow of ionic current arising from a transthreshold voltage step is due to the sodium ion influx. The activation process is characterized by a rapid increase in sodium permeability. The net driving force for sodium is the di↵erence between the transmembrane potential of 20 mV and the sodium Nernst potential of 57 mV, resulting in a driving force of 20–57 or –37 mV. Since this is negative, it is inward; consequently, a resultant inward (sodium) current is expected. This inward flow constitutes a bulge, because the elevated sodium permeability is transitory.


Figure 5.7: - Illustrative example of the ionic currents for a squid axon assuming the application of a voltage clamp of Vm = 20mV . Image Credits [2].

In fact, as the sodium permeability falls the potassium permeability rises and remains elevated. Assuming a potassium Nernst potential EK =-70mV, the potassium driving force is 20-[-70]=90mV and, since the driving force is positive, the current is outward. Valuable insight on the early membrane current can be achieved by clamping it to a value of Vm =EN a . For a short interval there is no current at all, even though in this early phase of the action potential, sodium permeability is tremendously higher. The reason is that there is no net force to cause a sodium current to flow. In Figure 5.8 the measured transmembrane currents arising from a series of voltage clamp experiments of di↵erent magnitude is shown. The figure includes a clamp at 117 mV (resting potential at -60mV, Vm =57mV), which corresponds to the sodium equilibrium condition, and we note the early measured current to be zero. The abolition of an early current when Vm is at the sodium Nernst potential confirms that it is the sodium ions that are responsible for this phase of total current. When the clamp voltage goes over 117 mV, the net driving force on sodium, (Vm

EN a ) is outward.

Note that for this condition the early current bulge is outward. The reversal potential is the value of voltage clamp for which the early inward current equals zero. From


the above argument, the reversal potential equals the sodium Nernst potential. This equality proves to be only a good approximation, as will now be explained.

Figure 5.8: - Measured Ion Currents for the squid axon following the application of a voltage clamp of the indicated value. The sodium Nernst potential is reached with a step change of 117mV because of a resting potential of -60mV and a Nernst potential of the sodium EN a =57mV. Image Credits [2].

NA/K ionic currents separation Once stated that is possible to record the total current which flows thru the membrane (which is the sum of boh the potassium, the sodium and the leakage current) the target is to separate the individual currents in order to develop a mathematical model for describe the behaviour of each ionic current. This separation was accomplished through the following procedure [2]. As shown in the Figure 5.8 the transmembrane current presents an early sodium current which is followed by a late, steady-state, outward current due mainly to potassium. A separation of the sodium and potassium currents is necessary in order to model the behaviour of each ionic component alone. The separation was accomplished with the preparation of a second extracellular medium by replacing 90% of the extracellular sodium by the inert element choline (the introduction of choline was simply to maintain isotonicity). Thereby, extracellular sodium concentration was reduced by a factor of 10. Initially, EN a =57mV; the reduction of extracellular sodium by a factor of ten should


lower the Nernst potential by 58log10 =58mV, so that the new value of EN a =-1mV. Accordingly, if the voltage-clamp experiment is repeated (with the same Vm = 56 mV) in the 10% sodium seawater, the sodium should be essentially in equilibrium. In this case the transmembrane current contains potassium only.

Figure 5.9: - Analysis of the Ionic Current in a Loligo axon during a voltage clamp. Trace A shows the response to a depolarization of 56 mV with the axon in seawater. Trace B is the response with the axon in a solution comprising 10% seawater and 90% isotonic chloride solution. Trace C is the di↵erence between traces A and B. Normal EN a =57mV, and in the reduced seawater EN a =-1mv. Image Credits [2].

Hodgkin and Huxley made a key assumption that the sodium and potassium ion fluxes are independent of each other (asserting the independence principle). In other words, they assumed that potassium current and sodium current crossed the membrane independently, each in its own pathways. Done this, Hodgkin and Huxley performed a series of voltage-clamp experiments for increasingly de-polarizing values. For each value of clamp voltage two experiments were performed. The first was for normal composition seawater and the second with a low-sodium seawater (replace 90% sodium chloride by choline chloride while potassium and remaining chloride ions are unchanged). In the second experiment, a potassium (only) current arises, so that, given the independence


principle, separation into IK and IN a is relatively easy. All the data collected by Hodgkin and Huxley from their voltage clamp experiments on the squid axon were the basis of a quantitative model for the squid axon membrane behaviour under both threshold and suprathreshold conditions which will results in the famous Hodgkin and Huxley mathematical model [17].

The Hodgkin and Huxley model Starting from experimental measurements of ionic currents, a mathematical model were constructed by Hodgkin and Huxley to fit all the data; their target was to write equations (and this was the most critical component of the model) to describe the single ionic currents in order to describe the results achieved with the voltage clamp. A remarkable aspect of the resulting equations is that they allow a membrane model to be constructed that not only reproduces the voltage clamp data itself but is capable of simulating a lot of phenomena, such as the generation and the propagation of action potentials. Hodgkin and Huxley considered the electric current flowing across the cell membrane during activation to be described by what we now call the parallel conductance model, which for the first time separated several ion-conducting branches. This model is illustrated in Figure 5.10. It consists of four current components: Current carried by sodium ions Current carried by potassium ions Current carried by other ions (designated leakage current, constituting mainly from chloride ions) Capacitive (displacement) current In this model, each of these four current components is assumed to utilize its own path or channel. To follow the modern sign notation, the positive direction of membrane current and Nernst voltage is chosen to be from inside to outside. The model is constructed by using the basic electric circuit components of voltage source, resistance, and capacitance as shown in Figure 5.10. The ion permeability of the membrane for sodium, potassium, and other ions is taken into account through the


Figure 5.10: - Electrical circuit representing membrane. RN a =1/gN a ; RK =1/gK ; Rl =1/gl . RN a and RK vary with time and membrane potential; the other components are constant. Image credits [17].

specification of a sodium, potassium, and leakage conductance per unit area (based on Ohm’s law) as follows: GN a =

GK =

GL =

IN a Vm VN a IK Vm







where: GN a , GK , GL = membrane conductance per unit area for sodium, potassium and other ions [S/cm2 ] IN a , IK , IL = the electric current carried by sodium, potassium and other ions [mA/cm2 ] VN a , VK , VL = Nernst voltage for sodium, potassium and other ions [mV ] Vm = membrane voltage [mV]


The conductance is a quantity that is determined experimentally by voltage clamp experiments. The major task of the HH model is to have developed equation which can describe correctly the ionic conductances without the need of a voltage clamp experimental session.

Potassium conductance: Because the behaviour of the potassium conductance during the voltage clamp experiment is simpler than that of the sodium conductance, it will be discussed first. The time course of the potassium conductance (GK) associated with a voltage clamp is described in Figure 5.11 [16] and is seen to be continuous and monotonic. Hodgkin and Huxley noted that this variation could be fitted by a first-order equation toward the end of the record, but required a third- or fourth-order equation in the beginning. This character is, in fact, demonstrated by its sigmoidal shape, which can be achieved by supposing GK to be proportional to the fourth power of a variable, which in turn satisfies a first-order equation. Hodgkin and Huxley gave this mathematical description a physical basis with the following assumptions. As is known, the potassium ions cross the membrane only through channels that are specific for potassium. Hodgkin and Huxley supposed that the opening and closing of these channels are controlled by electrically charged particles called n-particles. These may stay in a permissive (i.e., open) position (for instance inside the membrane) or in a nonpermissive (i.e., closed) position (for instance outside the membrane), and they move between these states (or positions) with first-order kinetics. The probability of an n-particle being in the open position is described by the parameter n, and in the closed position by (1 - n), where 0  n  1. Thus, when the membrane potential is changed, the changing distribution of the n-particles is described by the probability of n relaxing exponentially toward a new value.


Figure 5.11: - Behaviour of potassium conductance as a function of time in a voltage clamp experiment. The displacement of transmembrane voltage from the resting value [in mV] is shown (all are depolarizations). Image credits [16].

In mathematical form, the voltage- and time-dependent transitions of the n-particles between the open and closed positions are described by the changes in the parameter n with the voltage-dependent transfer rate coefficients ↵n and


This follows a first-

order reaction given by:

Figure 5.12: -

where ↵n = the transfer rate coefficient for n-particles from closed to open state [1/s] n

= the transfer rate coefficient for n-particles from open to closed state [1/s]

n = the fraction of n-particles in the open state 1 - n = the fraction of n-particles in the closed state


If the initial value of the probability n is known, subsequent values can be calculated by solving the di↵erential equation: dn = ↵n (1 dt




Thus, the rate of increase in the fraction of n-particles in the open state dn/dt depend on their fraction in the closed state (1

n) , and their fraction in the open

state n, and on the transfer rate coefficients ↵n and


Because the n-particles are

electrically charged, the transfer rate coefficients are voltage dependent. Furthermore Hodgkin and Huxley supposed that the potassium channel will be open only if four n-particles exist in the permissive position (inside the membrane) within a certain region. It is assumed that the probability of any one of the four n-particles being in the permissive position does not depend on the positions of the other three. Then the probability of the channel being open equals the joint probability of these four n-particles being at such a site and, hence, proportional to n4 . Finally, the potassium conductance per unit area is then the conductance of a single channel times the number of open channels. Its expression, if GKM ax is the conductance per unit area when all channels are open (i.e., its maximum value) is: GK = GKM ax n4


where GKM ax is the maximum value of potassium conductance [mS/cm2 ], and n obeys equation 5.8; The latter two equations, 5.8 and 6.10, are the basic expressions in the HH formulation.


Sodium conductance: The results obtained for sodium conductance in their voltage clamp experiments are shown in Figure 5.13 [16]

Figure 5.13: - Behaviour of sodium conductance as a function of time in a voltage clamp experiment. The displacement of transmembrane voltage from the resting value [in mV] is shown (all are depolarizations). Image credits [16].

The behaviour of sodium conductance is initially similar to that of potassium conductance, except that the speed of the conductance increase during depolarization is about 10 times faster. The rise in sodium conductance occurs well before the rise in potassium conductance becomes appreciable. Hodgkin and Huxley assumed again that at the sodium channels certain electrically charged particles called m-particles exist whose position control the opening of the channel. Thus they have two states, open (permissive) and closed (nonpermissive); the proportion m expresses the fraction of these particles in the open state (for instance inside the membrane) and (1-m) the fraction in the closed state (for instance outside the membrane), where 0  m  1. The

mathematical form for the voltage- and time-dependent transitions of the m-particles between the open and closed positions is similar to that for potassium. We identify


these with a subscript ”m”; thus the voltage-dependent transfer rate coefficients are ↵m and


These follow a first-order process given by

Figure 5.14: -

where ↵m = the transfer rate coefficient for m-particles from closed to open state [1/s] m

= the transfer rate coefficient for m-particles from open to closed state [1/s]

m = the fraction of m-particles in the open state 1 - m = the fraction of m-particles in the closed state An equation for the behaviour of sodium activation may be written in the same manner as for the potassium, namely that m satisfies a first-order process: dm = ↵m (1 dt




On the basis of the behaviour of the early part of the sodium conductance curve, Hodgkin and Huxley supposed that the sodium channel is open only if three m-particles are in the permissive position (inside the membrane). Then the probability of the channel being open equals the joint probability that three m-particles in the permissive position; hence the initial increase of sodium conductance is proportional to m3 . The main di↵erence between the behaviour of sodium and potassium conductance is that the rise in sodium conductance, produced by membrane depolarization, is not maintained. Hodgkin and Huxley described the falling conductance to result from an inactivation process and included it by introducing an inactivating h-particle. The parameter h represents the probability that an h-particle is in the non-inactivating (i.e., open) state - for instance, outside the membrane. Thus (1

h) represents the

number of the h-particles in the inactivating (i.e., closed) state - for instance, inside the membrane. The movement of these particles is also governed by first-order kinetics: where ↵h = the transfer rate coefficient for h-particles from closed to open state [1/s]


Figure 5.15: -


= the transfer rate coefficient for h-particles from open to closed state [1/s]

h = the fraction of h-particles in the NON-ACTIVATING state 1 - h = the fraction of m-particles in the INACTIVATING state This first order kinetic satisfies a similar equation to that obeyed by m and n: dh = ↵h (1 dt




Again, because the h-particles are electrically charged, the transfer rate coefficients ↵h and


are voltage-dependent and the sodium conductance is assumed to be pro-

portional to the number of sites inside the membrane that are occupied simultaneously by three activating m-particles and not blocked by an inactivating h-particle. Consequently, the behaviour of sodium conductance is proportional to m3 h, and GN a = GN aM ax m3 h


where GN aM ax is the maximum value of sodium conductance [mS/cm2 ], and m and h obeys, respectively, equation 6.11 and 5.11;

Transfer rate coefficients The transfer rate coefficients for ↵ and

for the gating variables n, m and h are deter-

mined from the next six equations which was developed by HH and, when substituted in equations 5.8, 6.10 and 5.11 lead to the curves plotted in figures 5.11 and 5.13. The dimension for the transfer rate coefficients is [1/ms]: ↵n =


0.1 e(1


0.01V 0 0.1V 0 ) 1 0.125

e(0.0125V 0 )




↵m =







e(V 0 /18) 0.07 e(0.05V 0 )

↵h =

In these equations V 0 = Vm

0.1V 0 )




0.1V 0



1 e(3

0.1V 0 )



Vr where Vr is the resting voltage. All voltages are given


in millivolts. Therefore V is the deviation of the membrane voltage from the resting voltage in millivolts, and its positive if the potential inside the membrane changes in the positive direction (relative to the outside). The equations holds for the giant axon of the squid at a temperature of 6.3 C. Note that in the voltage clamp experiment the ↵ and

are constants because the

membrane voltage is kept constant during the entire procedure. During an unclamped activation, where the transmembrane voltage is continually changing, the transfer rate coefficients will undergo change according to the above equations. Finally, is possible to summarize all the Hodgkin and Huxley equations [17]: Transmembrane Current and ionic pumps activation rates

I = CM

dVm + (Vm dt

VN a )GN a + (Vm

VK )GK + (Vm



dm = ↵m (1 dt




dn = ↵n (1 dt





dh = ↵h (1 dt




Ionic Conductances GN a = GN aM ax m3 h GK = GKM ax n4 GL = constant Transfer Rate Coefficients 0.01V 0 ↵n = e0.1 (1 0.1V 0 ) 1 2.5 0.1V 0 ↵m = e(2.5 0.1V 0 ) 1 0.07 ↵h = e(0.05V 0 )




0.125 0 e(0.0125V ) 4 m = e(V 0 /18) 1 0 e(3 0.1V ) +1


In addition to the variables discussed above, the constants of the Hodgkin-Huxley model are given here. The voltages are described in relation to the resting voltage (as shown): CM = 1 µF/cm2 EN a = V r EK = V r EL = V r

VN a =

115 mV

VK = +12 mV VL =

10.613 mV

GN aM ax = 120 mS/cm2 GKM ax = 36 mS/cm2 GL = 0.3 mS/cm2 Note that the value of VL is not measured experimentally, but is calculated so that the current is zero when the membrane voltage is equal to the resting voltage. These equation were solved numerically by Hodgkin and Huxley in order to validate them and to check if they can correctly predict the experimental evidences. The simulate a serie of voltage clamp experiments, at di↵erent depolarizations on the same axon as is possible to see in Figure 5.16; this validation step consists of ten depolarization on a giant squid axon at 4 C.


Figure 5.16: - Screenshot of the original HH article; the figure shows a comparison between calculated (analytically) and experimental Voltage Clamp data. Image credits [17].

Observing this last figure is possible to observe the goodness of the model and how the numerical data matches perfectly the experimental ones and this confirms that the solution of these equation can be useful for, as told before, simulate voltage clamp experiments without any practical experiments.


Extensions of the HH model One of the greatest advantages in the Hogkin and Huxley model is that, in spite it have been developed in order to describe the voltage clamp experimental observations, it can be applied - and ampliated - to successfully describe a lot of experimental setups and evidences, such as the generation of an action potential dued by an intracellular current pulse and the propagation of the signal along the axon. In the next subsections we will briefly see how the main HH equation can be modified to fully describe such situations.

Membrane action potential By a ’membrane’ action potential is meant one in which the membrane potential is uniform, at each instant, over the whole of the length of fibre considered [17]. There is no current along the axis cylinder and the net membrane current must therefore always be zero, except during the stimulus. If the stimulus is a short shock at t = 0, the form of the action potential should be given by solving equation 5.19 with I = 0 and the initial conditions that V = V0 and m, n and h have their resting steady state values, when t = 0. With appropriate instrumentation [16] is hence possible to stimulate an axon throughout its entire length with both a voltage or a current pulse; in this way the membrane voltage at each instant of time is identical over the entire length of the axon. This particular situation can be brought about by inserting a thin stimulation electrode along the axis of the entire fiber, whereas the other electrode, a concentric metal cylinder of the same length, is outside the axon. As a result, there is complete longitudinal uniformity of the potential along the axon. This means that the potential can vary only with respect to the radius from the axis, and only radial currents can arise. Furthermore, all membrane elements behave synchronously, so the entire axon membrane behaves as whole. Consequently, between the concentric electrodes, a membrane current will be measured that obeys the first equation in the HH model.


As for the voltage clamp, the two scientists Hodgkin and Huxley have made a comparison between calculated and experimental membrane action potentials; their results are reported in the figure 5.17:

Figure 5.17: - Screenshot of the original HH article; Upper traces: solution of the HH equation for initial depolarizations of 90, 15, 7 and 6 mv calculated at 6 C. Lower traces: membrane action potentials recorded at 6 C for giant squid axon. The numbers attached to the curves give the shock strength in mµcoulomb/cm2 . Image Credits [17]

The results here showed confirms the goodness of the HH model to simulate the generation of action potential as a consequence of di↵erent initial stimulus and its versatility.

Stimulation of excitable tissues The response of excitable cells to naturally occurring or artificial stimuli is a subject of great importance in understanding natural function of nerve and muscle, because most stimuli are produced by the natural system itself. Both electric and magnetic field stimulation are used in research investigations and in clinical diagnosis, therapy, and rehabilitation [2, 19].


Action potentials can be thus generated stimulating the excitable membrane with a current pulse injected at some time. In this experiment the user will measure the transmembrane voltage variation after an external proper current pulse. After a current pulse, the membrane voltage will not change immediately but will reach the final value in an exponential way (following the laws of the charging of a capacitor). With this technique is possible to stimulate the nervous cell and to generate an action potential (or, with long time pulses, a train of action potentials) and to check the time variation of the single ionic currents and the time variation of the gating variables m, n and h. In this experiment the neuron is impaled with one electrode to measure the resting potential and a second electrode called stimulating electrode. The stimulating electrode is connected through a switch to a battery. If the battery is oriented such that the positive pole is connected to the switch, closing the switch will make the inside of the cell somewhat more positive depending upon the size of the battery. (Such a decrease in the polarized state of a membrane is called a depolarization.) If we have the possibility to change the size of the battery (i.e. impress bigger current intensity), initially, the switch closure produces only small depolarizations. However, increasing the size of the battery, the potentials become larger and eventually the depolarization is sufficiently large to trigger an action potential. The action potential is associated with a very rapid depolarization to achieve a peak value of about +40 mV in just 0.5 milliseconds (msec). The peak is followed by an equally rapid repolarization phase. This kind of experiment is very useful because it permits to observe the typical everything-or-nothing behaviour of the excitables membranes applying under and over threshold pulses. The voltage at which the depolarization becomes sufficient to trigger an action potential is called the threshold. If a larger battery is used to generate a suprathreshold depolarization, a single action potential is still generated and the amplitude of that action potential is the same as the action potential trigged by a just-threshold stimulus.


Figure 5.18: - Di↵erent membrane responses at di↵erent current stimulus. Image Credits [19].

The sample recording in the latter figure illustrates two very important features of action potentials. First, they are elicited in an all-or-nothing fashion. Either an action potential is elicited with stimuli at or above threshold, or an action potential is not elicited. Second, action potentials are very brief events of only about several milliseconds in duration. Initiating an action potential is somewhat analogous to applying match to a fuse. A certain temperature is needed to ignite the fuse (i.e., the fuse has a threshold). A match that generates a greater amount of heat than the threshold temperature will not cause the fuse to burn any brighter or faster. Just as action potentials are elicited in an all-or-nothing fashion, they are also propagated in an all-or-nothing fashion. Once an action potential is initiated in one region of a neuron such as the cell body, that action potential will propagate along the axon (like a burning fuse) and ultimately invade the synapse where it can initiate the process of synaptic transmission. In the example in the figure, only a single action potential was generated because the duration of each of the two suprathreshold stimuli was so brief that sufficient time was only available to initiate a single action potential (i.e., the stimulus ended before the action potential completed its depolarization-repolarization cycle). But, as shown in next figure, longer-duration stimuli can lead to the initiation of multiple action potentials, the frequency of which is dependent on the intensity of the stimulus. Therefore, it is evident that the nervous system encodes information not in terms of the changes


in the amplitude of action potentials, but rather in their frequency. This is a very universal property. The greater the intensity of a mechanical stimulus to a touch receptor, the greater the number of action potentials; the greater the amount of stretch to a muscle stretch receptor, the greater the number of action potentials; the greater the intensity of a light, the greater the number of action potentials that is transmitted to the central nervous system. Similarly, in the motor system, the greater the number of action potentials in a motor neuron, the greater will be the contraction of the muscle that receives a synaptic connection from that motor neuron. Engineers call this type of information coding pulse frequency modulation. Hodgkin and Huxley did not validate this kind of experiment but is possible to find works on literature in which this has been done [20]. To simulate this kind of experiments the main equation of the HH model has to modified in order to simulate the presence of the current pulse: CM

dVm = dt

Iionic + Ist


in which the term Ist represent the current pulse expressed in µA/cm2 and it can be a time varying function such as a step function; Except this one the other equations of the HH model remains the same. The next step is to check how this mathematical model can be slightly modified in order to predict not only the generation but also the propagation of action potential along an active fiber.

Propagating nerve pulse Finally, to conclude this chapter entirely dedicated on the Hodgkin and Huxley model, let see how it can successfully simulate the more complex and near to reality situation; the propagation af o nervous pulse. If a long thin fiber is initially depolarized at one end, the initial depolarization at first has a limited extent. Thereafter propagation from the active (already depolarized) region to adjoining regions will take place. In this context the term propagation


means that each patch of excitable membrane initiates an action potential at an adjacent patch, which then creates its own sequence of action potential events, including initiating yet another patch. This process is quite complex in that it involves multiple things happening at the same time. In the propagated action potential the local circuit currents have to be provided by the net membrane current; with the aid of the core-conductor model [2] and many other approximations, this lead to the relation: i=

1 @ 2 Vm r1 + r2 @x2


in which i is the membrane current per unit length, r1 and r2 are the extracellular and intracellular resistances (per unit length), and x is the distance along the fibre. If the axon is surrounded by a large volume of conducting fluid, r1 will be negligible with respect to r2 with the result i=

1 @ 2 Vm r2 @x2


a @ 2 Vm 2R2 @x2


or I=

in which I is the membrane current density, a is the radius of the fibre and R2 is the specific resistance of the axoplasm. If now we insert this last expression in the main HH equation the result will be: @ 2 Vm @Vm = Cm + (Vm 2 @x @t

VN a )GN a + (Vm

VK )GK + (Vm



with all the subsidiary equation unchanged [2]. The membrane current equation 6.4 relates the density of current crossing the membrane (i.e. the transmembrane current) to the second spatial derivative of the transmembrane potential. From a physiological perspective this equation is the result of that in a fiber the current can flow axially. Regions are linked by currents coming inward across the membrane at one place, and moving outward at another. In other words, the mechanism that allows an action potential at one place to trigger an action potential at a site nearby is the mechanism of axial currents from one place to another. The result is that current exits the membrane at places di↵erent from where it enters,


leading to excitation at a di↵erent membrane location. As the original HH equations, the relation 6.4 is a partial di↵erential equation which cannot be solved analytically and requires particular numerical methods to solve it as we will see in the next chapter in which we will validate our numerical environment in order to reproduce the original results obtained by Hodgkin and Huxley and then we will apply all our know-how to extract important information about the relation between intracellular and extracellular potentials near the active fibre.

Inclusion of the temperature e↵ects Adaptations to changes in the external temperature are an essential component of living systems, as most physiological systems function at di↵erent rates at di↵erent temperatures. The expressions given by Hodgkin and Huxley for the rates ↵ and rate constants are for a temperature of 6.3 C, which is a natural temperature for the squid (remember that the HH model have been derived from studies on the axon of the giant squid) in seawater. Rate constants vary with the temperature and that changes are taken in account by adjusting the gating variables expressions (which are function of the ↵ and

values) as follow [20]:

8 dn > > ( = ↵n (1 n) n n) · k > > > < dt dm (5.28) ( = ↵m (1 m) m m) · k > dt > > > > : ( dh = ↵ (1 h) h h h) · k dt The coefficient k here inserted in the equations for the gating variables is called thermic coefficient and takes in account the e↵ect of the temperature and is expressed as: (T T0 )/10

k = Q10


in which T is the actual temperature, T0 is the temperature in which the HH model was derived, 6.3 C, and the parameter Q10 was set equal to 3 by Hodgkin and Huxley as a good value for simulate correctly the behaviour of the squid membranes at di↵erent temperatures (but not over 31 C) [20, 21].


References [1] Berne, R.M., Levy, M.N. Levy, Fisiologia, V Edizione, Casa Editrice Ambrosiana. [2] Plonsey, R., Barr, R., Bioelectricity: a quantitative approach, (2007), Third Edition, Springer. [3] Hille, B. Ionic Channel of Excitable Membranes, (1984), Sinauer Associates, Sunderland, MA, USA; 1st Edition. [4] Bullock T. H., Orkand R., Grinnell A., Introduction to Nervous System, (1977), W. H. Freeman ed. [5] Purves D., Augustine G. J., Fitzpatrick D., Katz L. C., LaMantia A-S., McNamara J. O., and Williams S. M., Neuroscience, 2nd edition, (2001), Sunderland (MA): Sinauer Associates. [6] Junge D., Nerve and Muscle Excitation, 2nd edition, (1981), Sunderland, MA: Sinauer Associates. [7] MacDonald P. E., Rorsman P., Oscillations, intercellular coupling, and insulin secretion in pancreatic beta cells, (2006), PLoS. Biol., 4 (2): e49. [8] http : //hyperphysics.phy

[9] http : // otential [10] Young JZ (1936), The giant nerve fibers and epistellar body of cephalopods. Q. J. Microsc. Sci. 78: 367-86. [11] Cole, K.S., Dynamic Electrical characteristics of squid axon membrane, (1949), Arch. Sci. Physiol. 3:253-8.


[12] Marmont G., Studies on the axon membrane. A new method., (1949), J: Cell. Comp. Physiol. 34:351-82. [13] Sakmann B, Neher E (1984), Patch clamp techniques for studying ionic channels in excitable membranes. Annu. Rev. Physiol. 46: 455-72. [14] http : // rizes/medicine/laureates/1991/ [15] Hodgkin A. L., Huxley A. F., Katz B., Measurement of current-voltage relations in the membrane of the giant axon of Loligo., (1952), J. Physiol. (Lond.), 116: 424-48. [16] http : //www.bem.f i/book/04/04.htm [17] Hodgkin, A. L., Huxley, A. F., A quantitative description of membrane current and its application to conduction and excitation in nerve., (1952), The Journal of Physiology, 500–544. [18] Hodgkin A., A note on Conduction Velocity, (1954), J. Physiol. (Lond.), 125:221-4. [19] http : //, John H. Byrne, Ph.D., Department of Neurobiology and Anatomy, The UT Medical School at Houston. [20] Rattay, F., Aberham, M., Modeling axon membranes for functional electrical stimulation, (1993), IEEE Transactions on Bio-Medical Engineering, 40(12), 1201–9. [21] Rattay F., Electrical Nerve Stimulation, (1990), New York: Springer.



FEM simulation of Action Potentials In the previous chapter we have seen how is possible to reduce the generation and propagation of nervous signals to a PDEs based mathematical model known as HH model; in this chapter we will show how we have integrated this model into a numerical environment (COMSOL Multiphysics software) in order to solve and validate it as already done for the Fokker-Planck equation in the case of the dynamic a rotary molecular motor (Chapters 3 and 4). This validation step will start from the original article in which the mathematical model is presented and compared to experimental results in order to be sure that our numerical results agrees them. Done this, we will exploit the COMSOL capacities to couple model governed by di↵erent physics in order to compute the extracellular potentials near an active nervous fibre; we will hence focus our attention on extracellular action potentials.

Validation of our FEM environment Validate a numerical model means to be sure that the results are in good agreement with experimental evidences; we will hence solve the HH model for the basic experimental setups like the voltage clamp and the space clamp and we will check the results with the one obtained in original HH work.


The first thing to do is to insert the HH mathematical model into our numerical environment; to simulate the Hodgkin and Huxley model the COMSOL Multiphysics software have been used (see chapter 2 for a description and reference of this software). The four di↵erential equations of the HH model have been inserted in COMSOL modifying the pre-existent coefficient form PDE which is already inserted in the program as done for the FPE equation; done this, a proper domain has to be prepared. For this purpose a simple, monodimensional, linear domain have been created. The equations will be solved on this line which represents the axon on which the experiments are done. At the two side of this line the default Neuman boundary conditions have been set in order to allow the mathematical function to have a real, non-zero, value on these boundaries.

Numerical simulations of the Voltage Clamp First step of this validation stage is to check how the numerical results matches the experimental voltage clamp results [1]. In previous chapter we have seen that in a voltage clamp experiment the membrane voltage is keeped constant along the nerve during all the experiment; mathematically, we have: dVm =0 dt At a constant voltage the coefficients ↵n ,


(6.1) ↵m ,


↵h and


are constant and

the solution is obtained directly in terms of the equation already given for m, n and h. The total ionic currents was computed from these for a number of di↵erent clamp voltages and is compared with a series of experimental data and analytic solutions of the HH system. With in mind the figure 5.16 of the previous chapter, the results obtained computing the HH equations in COMSOL with the same parameters are shown in Figures 6.1, 6.2, 6.3 and 6.4. Comparing the result obtained by our numerical simulations with the experimental ones [2], confirms the goodness of the solutions obtainable by our FEM environment and let us go over with this validation step and to simulate the generation of action potentials.


Figure 6.1: - Computed membrane current after two voltage clamp depolarization of 5mV and 10mV.

Figure 6.2: - Computed membrane current after two voltage clamp depolarization of 10mV and 30mV.


Figure 6.3: - Computed membrane current after two voltage clamp depolarization of 40mV and 60mV.

Figure 6.4: - Computed membrane current after two voltage clamp depolarization of 80mV, 100mV, 115mV and 130mV.


Membrane Action Potential Generation Once validated the model for the Voltage Clamp is time to go over and to check how is possible to simulate the generation of an action potential along a fiber; such simulation have an analogue experimental setup in the space clamp technique [3]. This technique allows the stimulation of an axon throughout its entire length; making the membrane voltage, at each time instant, identical over it; as a result, there is no current along the axis cylinder and the net membrane current must therefore always be zero, except during the stimulus. Consequently, if the stimulus is a short shock at t = 0, the form of the action potential should be given solving the HH model with Im = 0; Cm dVm /dt =

II and the

initial conditions that V = V0 and m, n and h have their resting steady state values when t = 0. Validation of our numerical environment proceed solving the HH model in order to check if it can predict the generation of action potential in a way which is in perfect accordance with the experimental values and with the analytical solution of the HH model. We have thus solved the HH model following the experimental data; the main equation of the model which have been solved is: Im = 0,


dVm = dt



The other equations and parameters will be the same as the one used in the previous section. Di↵erent initial depolarizations have been tested in order to simulate the experimental results; four di↵erent simulation have been run with, as done in the original HH article, initial depolarizations of 90, 15, 7 and 6mV for a duration of 6ms [2].


Figure 6.5: - Computed membrane potentials on the giant squid axon for four di↵erent depolarizations. Is interesting to note that is possible to simulate the threshold mechanism; passing from 6mV to 7mV we have the activation of the fiber.

Solving numerically the HH model allows the user to obtain not only the time course of the membrane voltage but also the time course of every variable present in the model. Of particular interest are the single ionic currents and the ionic gating variables. The time course of the variables m, n and h is obtainable directly plotting, from Comsol, the correspondent variables; they will result in the plot of figure 6.6 which confirms the typical behaviour of the ionic pumps when an action potential is elicited as explained in the second chapter [4, 5]. As for the time course of the gating variables is possible to plot the time course of each ionic current as showed in figure 6.7. Finally, the sum of the curves represented in 6.7 is the total transmembrane current, showed in Figure 6.8.


Figure 6.6: - Time course of the three gating variables m, n and h which represents, respectively, sodium, potassium and leakage ionic pumps. These curves referee on a 7mV space clamp depolarization experiment.

Figure 6.7: - Time course of the three ionic components of the total transmembrane current. Also in this case the curves referee on a 7mV space clamp experiment.


Figure 6.8: - Time course of the total transmembrane current for a 7mV space clamp experiment.

This last quantity represents the total current which flows through the membrane; positive current is from inside to outside, negative represents the opposite. This quantity will be relevant in the next steps of this study because this quantity is the one which will be fundamental in relating the transmembrane potential to the extracellular potential caused by the nervous activity of the active fibres [4].

All these results, as for the voltage clamp validation steps, confirms another time the goodness of the model [1, 4, 5]. In the next section we will simulate the stimulation of an active fiber with di↵erent current pulses in order to show the possibility to simulate the typical everything-or-nothing behaviour of the active fibres.

Stimulation of active fibres with current pulses. In the previous chapter we have seen how is possible to extent the HH model in order to simulate the electric stimulation of active fibres; to reach this target we have seen that the equation for the transmembrane potential can be modified, including the


stimulating pulse: CM

dVm = dt

Iionic + Ist


in which the term Ist represent the current pulse expressed in µA/cm2 and it can be inserted as a function which vary in time such as a step function; Except this one the other equations of the HH model remains the same [4, 5]. This kind of simulation is very useful as it permits to observe the typical everythingor-nothing behaviour of the excitable membranes applying under -and over - threshold pulses [6]. In this validation step the current pulse will be a train of identical long pulses with di↵erent intensity in order to show how the system behaves in both the under- and over- threshold situations and in order to determine which is the threshold for this particular kind of axon (another time we are simulating the axon of the giant squid). In the figure 6.9 this behaviour is showed. For such an experiment is possible to analyse the relation between the strength of the current stimulus and its duration in relation to his ability to generate an action potential. For the HH model of the giant squid axon a number of experiments have been carried out with the purpose of know how long and intense a current stimulus should be in order to generate an action potential. In table 6.1 we show the strength-duration relation for five di↵erent cases; for every line of the table, a current clamp experiment with di↵erent stimulus amplitude and duration has been carried out. Table 6.1: Strength-Duration relationship between current pulses and the duration in order to elicit an action potential. Duration [ms]

Strength [µA/cm2 ]












Figure 6.9: - Time course of the membrane potential after three - equal length - current stimulus pulses at increasing amplitude. The one which elicits an AP is at 72 µA/cm2 for 100 µ s.

Its also interesting to note that is possible to simulate the refractory period; after an action potential is elicited, if a second pulse is injected while the system is in this phase, the membrane will not be re-excited during the so-called refractory state. This behaviour is showed in the figure 6.10 as a result of an opportunely designed simulation. Plotting the time course of the gating variables related to this simulation can be of great help understanding refractoriness; observing the plots in figure 6.11 one can appreciate that refractory period is manly governed by the behaviour of the inactivating parameter h. Following an action potential, h decreases to a very low value. Consequently, even a very large stimulus elicits only a small sodium current, preventing re-excitation of the membrane.


Figure 6.10: - Time course of the transmembrane voltage after a train of two - equals in length and strength - current stimulus pulses. The second is injected while the system rely in his refractory phase and does not generates AP.

Figure 6.11: - Time course of the gating variables m, n and h after a train of two equals in length and strength - current stimulus pulses. The second is injected while the system rely in his refractory phase. Is possible to observe the particular behaviour of the h parameter which leads the total refractory period.


Is also interesting to note that the HH model, which, as said, is derived for the giant squid axon, can be modified in order to completely simulate the behaviour of many di↵erent classes of excitable cells of many di↵erent species [7]. As for an example we report here the simulations of the generation of action potentials for the axon of the rat subthalamic nucleus neuron. The results obtained with such a simulation are obtained simulating a current-clamp experiment injecting a current pulse of 10 µA/cm2 in the middle of a 9.4 µm radius fibre for a duration of 100ms; the time course of the generated action potential is showed in figure 6.12 ad the results are in good agreement with experimental evidences [7, 8].

Figure 6.12: - Train of action potentials generated on the rat subthalamic axon after a 10 µA/cm2 current pulse for 100 ms

Another time, our numerical environment, in combination with the HH model demonstrate its power and its great usefulness in simulating real cells behaviours and understanding them. In order to complete this first part of this chapter, entirely dedicated to the numerical simulations of the HH models, let see the propagation of action potentials.


Propagating nerve pulse Starting from the last subsection is possible to analyse how an action potential, after its generation with an opportune over threshold current pulse, propagates along the fibre axis. In the previous chapter we have seen that the HH model can be modified to successfully simulate the signal propagation; with such a modification the equation 6.4 describes the space and time variation of the transmembrane potential and all the subsidiary equation being unchanged.

@ 2 Vm @Vm = Cm + (Vm @x2 @t

VN a )GN a + (Vm

VK )GK + (Vm



The action potential is generated in the point where the current pulse is injected; after this, it will propagate along the fibre. In the next figure the space variation of the membrane potential is shown; the plot contains two curves, one for the membrane potential at time t = 4ms and the second at time t = 5ms. From this plot is possible to obtain the propagation speed (in our simulations the obtained propagation velocity is equal to 19 m/s) which is very similar to what was obtained by Hodgkin and Huxley in their original work [2, 9], 18.8 m/s.

Figure 6.13: - Space course of the membrane potential at the times t = 4ms and t = 5ms.


In the last section we have seen how is possible to modify the HH method in order to describe the electric behaviour of fibres for di↵erent kind of species; in the same way is possible to simulate the propagation of the signals along such a fibres. As before we have simulate the stimulation of a rat subthalamic axon injecting a current pulse of 10 µA/cm2 in the middle of a 9.4 µm radius fibre for a duration of 100ms; plotting two consecutive action potentials allows us to calculate a propagation velocity of 0.5 m/s, which is very similar to experimental values. [8]

Figure 6.14: - Space variation of the transmembrane potential on the rat subthalamic axon after a 10 µA/cm2 current pulse for 100 ms for two consecutive time, 64 ms and 65 ms. From this plot one can extract the propagation velocity of 0,5 m/s

This last test of the HH model in combination with the Comsol FEM environment allows to a↵ord that the model can fully simulate the typical behaviours of an excitable tissue; it can successfully simulate the generation of an action potential and its propagation along the fibre not only for a giant squid axon (the HH model have been designed to describe this nerve) but also we can simulate axon and cells of di↵erent species such as the rat subthalamic axon. In spite of all these results, one of the limits of the HH


model is its monodimensionality; it can successfully model only what happen on the fibre membrane and it doesn’t tell us nothing on what happen in the space around it. In order to achieve more information about this point, i.e. compute the extracellular electric potentials related to an action potential, we have to couple the HH model with an opportunely created model which can simulate the extracellular environment.

Computing Extracellular action potentials The extracellular field potential is the electrical potential produced by cells, e.g. nerve or muscle cells, outside of the cell; extracellular action potential (EAP) recordings form one of the primary means for study of the nervous cells behaviour [10, 11]. Simulation of extracellular fields is one of the substantial methods used in the area of computational neuroscience. Our main purpose here is to provide a computational model that can compute the spatio temporal variations of the extracellular electric fields produced by the action potentials generated by a current-clamp experiment on an axon. The extracellular potential induced by a spike in a neuron was calculated in two distinct stages. Firstly the Hodgkin and Huxley (HH ) model is used for simulate the generation and propagation of the action potential along an axon (squid axon) then the membrane currents will be used to compute the extracellular potentials as described below.

Theoretical background From the first principle of the Maxwell’s theory of the electromagnetism, the electric properties of a conductive medium are fully determined by two parameters, conductivity

and permittivity ✏ [12]; although conductivity quantifies the local relation between

the electric field and the current, permittivity characterizes the response of the system in terms of separation of opposite charges (polarization) in the presence of an electric field. Maxwell’s theory of electromagnetism allows one to compute electric and magnetic fields or potentials, for a given distribution of charges and electric currents. Because charges move very slow in biological media, the e↵ects of magnetic fields are


very small compared to that of the electric field and will be neglected here. One of Maxwell’s equations is Gauss law: r · (✏E) = ⇢


Here E denotes the electric field, D=✏E is the displacement, and ⇢ is the charge density of the extracellular medium, also allowed to vary slowly in time. The continuity equation relates the current density j to the charge density ⇢: @⇢ =0 (6.6) @t which states a balance between the electric flux into some volume and the change r·j+

of the total charge in this volume. In other words, no charge will get lost (continuity). Finally, combining Ohms law j = ( E)


with the continuity equation we have, for scalar conductivity: @⇢ =0 (6.8) @t In the following, we will assume that electric currents are distributed on the surface r · ( E) +

of the membrane (ionic currents), and that these currents are allowed to vary in time. Also, we have E =

rV because the magnetic field is negligible.

Model coupling Let see now how we have connected the action potential which is generated and propagates along the fibre to the extracellular field. A cylindrical fibre carrying an action potential in the x direction produces potentials throughout the surrounding medium. A sketch of the geometry used in the simulations is in the figure 6.15


Figure 6.15: - An action potential is propagating on a fibre in the positive x direction. The fibre is divided into mathematical segments (we’re using the FEM method for solving the equations), and one such segment is drawn. The fibre lies in a uniform extracellular medium of conductivity e . The field arising from the action currents at an arbitrary field point P is desired. Image Credits [2]

The fibre is divided conceptually into small elements along its length. Fiber element dx, identified in figure 6.15, lies within the region occupied by an action potential. Out of this di↵erential fibre element a current emerges into the extracellular region. The amplitude of this current from the element is the transmembrane current per unit length, im , times the length dx, i.e. im dx. For simplicity we have assumed the fibre to be thin relative to the extracellular volume, so that it can be treated mathematically as a line.

In general, for a current point source, the potential is given by:



1 4⇡


I0 r

where r is the distance from the point source to the field point,

(6.9) e

is the extracellular

electric conductivity, and I0 is the source current amplitude. In the case of a cylindrical fibre, the analogue to the point source I0 is the transmembrane current Im . The transmembrane current is, of course, not confined to a point, but instead varies along


the fibre. Therefore the contribution of im dx to


can be rewritten, according to the

last eq.: d



1 4⇡

im dx e r


Finding now the potential requires an integration over the full length of the fibre: 1 e (P ) = 4⇡




im dx r


Is important to note that the expression im dx is used instead of I0 , since im dx is the current emerging from one segment of the fibre (as shown in fig 6.15). The equation 6.10 is a di↵erential contribution to the extracellular potential d


To get


an integral over the length of the fibre is required, so that all segments along the fibre that have non-zero membrane current are taken into account. From the equation for im is possible to determine the final form for the extracellular potential


The membrane current, obtainable by the so-called core-conductor model,

inserted in the eq. 6.11 leads us to:


a2 = 4

i e


@Vm /@x2 dx r


The last equation is used widely as it depends on Vm , which often is known. This definition may be used to find potentials at field points that may be chosen to be either close to or far from the active fibre. The same equation can be used to find the potential at only one position, or repeatedly to find potentials for each location within a family of positions.

Computation of the extracellular potentials. results.

Model validation and

In the simulation environment the extracellular space is represented as a monodimensional square with homogeneous electric conductivity

and permittivity ✏. The nervous

fibre is represented as a line (on which the Hodgkin and Huxley model is valid and so solved) that acts as a current source (i.e. the transmembrane current). In other words this model needs a coupling between the electrodynamic problem of the extracellular


environment and the neurophysiologic problem of what is occurring on the nervous membrane. Our study on the extracellular potentials is carried out for the simplest case i.e. homogeneous extracellular medium characterised by a relative electric permittivity ✏r = 80 which is the relative permittivity of the water at 20 C and the electric conductivity is set equal to


= 1/330⌦cm in order to have results comparable with

other results present in literature [2, 12, 13, 14]. In the simulation environment, the protocol is to first solve the HH equation; from these extract the transmembrane currents and use them in the bidimensional model to evaluate the potentials. The extracellular action potentials computed in this section referee to the generation of an action potential on a giant squid axon after a current pulse of 90 µA/cm2 for a duration of 1ms. In the figure 6.16 the action potential generated by this pulse is shown. In this plot the time variation of the membrane potential at the point x = 4 mm is shown and the current pulse is injected at the point x = 1mm. This kind of simulation can also be carried out for every kind of cell just with a proper modification of the relative HH model.

Figure 6.16: - Time course of the membrane potential at the point x = 4mm.


Figure 6.17: - Total transmembrane current related to the action potential of fig. 6.16

In figure 6.17 the total transmembrane current related to the action potential of fig. 6.16 is plotted. This current is the one which will be coupled to the bidimensional domain and which is the responsible of the extracellular electric field generated by the action potential. Once the HH model is solved and the total transmembrane current is extracted is thus possible to go over with the bidimensional domain; sample output from the simulation environment is showed in the figure 6.18 where the total extracellular potential arising from the membrane at the time t=3.3ms is showed. This time moment is, approximately, the time at which the extracellular potential rise its maximum value for the point x=4mm. The line which appears in the middle of the figure is the monodimensional domain related to the fibre on which the HH model is valid and which acts as a current source density.


Figure 6.18: - Extracellular potential, in µV arising from the action potential which is running on the fibre membrane represented as a line.

From this first plot is possible to appreciate the intensity of the extracellular fields near the active fibre; at the action potential maximum we obtain an extracellular potential of about 109 µV while at its minimum (where the membrane potential is negative) we have a negative value of about -25µV. These values are in accordance with other experiments [12]. More detailed informations on the extracellular potential form and intensity can be extracted plotting the time course of the extracellular potential at di↵erent distances from the active fibre axis; this is what is shown in figure 6.19.


Figure 6.19: - Time course of the extracellular potential, in µV arising from fibre at di↵erent points which are at di↵erent distances from the fibre itself.

From this last plot is firstly possible to notice that, when the membrane potential is at its maximum (⇡ 40mV ) the extracellular is ⇡ 0.1mV and this value falls down in the extracellular space as more as we go far from the fibre; the extracellular potential has an average dimension order of the microvolts in great accordance with experimental recordings [12, 13, 14]. In the space distribution of the potential plays a central role the electrical conductivity of the medium; a higher value of this parameter will results, obviously, in a higher value of the EAP far from the fibre. Here we have used a value of this parameter in order to, as told before, simulate a typical physiologic extracellular environment.

From the simulation is finally possible to plot the space variation of the EAP at di↵erent distances from the active fibre. The results showed in figure 6.20, which refers to the time step 3.3ms, shows the typical waveform of the extracellular action potential


[13] and a maximum value of 80 µV , also in great accordance with experimental values [12, 13, 14].

Figure 6.20: - Space course of the extracellular potential, in µV arising from fibre at di↵erent distances from the fibre itself.

Results showed here demonstrates that with our approach is possible to successfully compute not only solutions of the HH models, but also we can couple them with the equations for the electric potential and obtain the extracellular potential generated by the membrane action potential. This result is of extreme importance; in spite of all it great features, the HH model has the problem that it cannot go over the membrane. It can successfully simulate all that happens on the membrane like transmittance potential and current and their variations both in space in time but it doesn’t tell us nothing of what happen around the active cell/fibre. The model here showed solves this problem allowing us to compute the extracellular potentials generated by an action potential


and to extract their form, intensity and variability both in space and time. Our model also guarantee us a great freedom; we can change the HH model associated to the fibre simulating di↵erent kind of fibres and cells and we can simply modify the parameters related to the extracellular environment in order to simulate how they changes the distribution of the extracellular potentials.

Conclusions In this chapter we have seen how a model created to simulate the generation and propagation of a nervous signal along an active fibre works great matching perfectly the experimental results obtained by Hodgkin and Huxley. This first result is of great importance because solutions of the HH model have never been obtained in such a simple way. There are software appositely designed to simulate the generation and propagation of the nervous signals using the HH model [15] but fully understand how they work, how to set up simulations and how to modify them in order to simulate di↵erent cells is quite complicated, also their interface is not user friendly; with our engineering approach we have set up a simple numerical environment, with an easy interface, which can give to the user solutions to the HH model in a few minutes with not a huge neurophysiological know-how. This model have been validated and all our results agrees the experimental evidences and the analytical solutions of the HH model. In the second section of this chapter we have exploited all the qualities of COMSOL Multipysics coupling a model described by the HH model with one described by the Maxwell’s equations for the electric potential in order to compute the extracellular action potentials. This result allows us to overcome the monodimensionality of the HH model and it can successfully describe what happen in the space around the active cell/axon membrane. The model described here is extremely simple and can be easily ampliated in order to simulate and describe all the side-e↵ects related to the extracellular action potentials like, one of the most important of them, the ephaptic coupling between adjacent active fibres [12, 16]. Ephaptic coupling is a form of communication within the nervous system and is distinct from direct communication systems like electrical synapses and chemical synapses. It can refer to the coupling of adjacent (touching) nerve fibres as a result of local electric fields. Ephaptic coupling


can influence the synchronization and timing of action potential firing in neurons [17]. Modelling this phenomena can be one of the first ampliation and application of this model; we are just working on a model which can simulate the extracellular stimulation of active fibres [4, 5]. Once completed, this model will be coupled with the one for the computing of the extracellular action potentials in which one fibre will act as a current source and the other will be stimulated by the electric potential produced by the first. In such an experiment great importance will be covered by the electric properties of the extracellular media [12].

References [1] Cole, K.S., Dynamic Electrical characteristics of squid axon membrane, (1949), Arch. Sci. Physiol. 3:253-8. [2] Hodgkin, A. L., Huxley, A. F., A quantitative description of membrane current and its application to conduction and excitation in nerve., (1952), The Journal of Physiology, 500–544. [3] Marmont, G. (1949). Studies on the axon membrane: I. A new method J. Cell. and Comp. Physiol. 34:351-382. [4] Plonsey, R., Barr, R., Bioelectricity: a quantitative approach, (2007), Third Edition, Springer. [5] Rattay, F., Aberham, M., Modeling axon membranes for functional electrical stimulation, (1993), IEEE Transactions on Bio-Medical Engineering, 40(12), 1201–9. [6] http : //, John H. Byrne, Ph.D., Department of Neurobiology and Anatomy, The UT Medical School at Houston. [7] Pospischil M., Rodriguez M.T. et al., Biol Cybern (2008) 99:427-441


[8] Stuart G., Schiller J., Sakmann B., Action potential initiation and propagation in rat neocortical pyramidal neurons, (1997), The Journal of Physiology, 505, 617-632. [9] Hodgkin A., A note on Conduction Velocity, (1954), J. Physiol. (Lond.), 125:2214. [10] Pettersen K. H., Lind´en H., Dale A. M., and Einevoll G. T., Extracellular spikes and current-source density, Handbook of Neural Activity Measurement, edited by Romain Brette and Alain Destexhe. [11] Nunez P. L., Srinivasan R., Electric Fields of the Brain., (2006), Oxford University Press. [12] Holt G. R., Koch C., Electrical interactions via the extracellular potential near cell bodies., (1999). Journal of Computational Neuroscience, 6(2), 169–84. [13] Gold C., Henze D. A., Koch C., Buzs´aki G., On the origin of the extracellular action potential waveform: A modeling study., (2006)., Journal of Neurophysiology, 95(5), 3113–28. [14] Pettersen K. H., Dale A. M., Einevoll G. T., Extracellular spikes and currentsource density, (2010), Handbook of Neural Activity Measurement, edited by Romain Brette and Alain Destexhe. [15] http : // [16] Binczak, S., Eilbeck, J., Scott, A. (2001). Ephaptic coupling of myelinated nerve fibers, (2001), Physica D: Nonlinear Phenomena, 1–26. [17] http : // oupling



Discussion and Future Perspectives In this Ph.D work we have focused our attention on two extremely di↵erent chemical and physical processes but both representable by opportunely designed mathematical models. To obtain important informations about these processes these mathematical models have to be solved. Since these models are based on partial di↵erential equations (Chapter 2) solution cannot be achieved analytically and we need a proper numerical environment in order to solve them. For this purpose we used the COMSOL Multiphysics software which uses the Finite Element Method (FEM) to get solutions to problems described by PDEs. This particular software was born for engineering applications and one of our main target was to adapt it in order to simulate systems of molecular dimensions. In the first section of this work we have seen how the dynamic of a rotary molecular motor can be successfully described by the Fokker Planck equation and solutions obtainable by it (Chapter 3) [1]. Firstly we used a traditional computational chemistry approach to obtain the energetic path in which the dynamic of such a molecule is exploited; this stage is fundamental as the Fokker-Planck equation has to be defined inside it; done this, the equation have been solved and, after a numerical elaboration of the results from it, we have extracted a lot of useful and interesting datas (Chapter 4).


First of all the determination of the energetic path gave us results in perfect accordance with the experimental evidences obtained by transient absorption spectroscopy analysis. Solving then the FPE, allowed us to fully understand the thermodynamic behind the process obtaining plots for the time course of both the internal energy, the entropy and the free energy changes during the process. Using the statistical thermodynamic tools to combine that results, we achieved a complete description of the dynamic of the molecular motor, obtaining at the end the total force exerted by the molecule during the transformation with results which are agree with the experimental evidences. We were also able to extract precious informations such as the total work done by the motor and the efficiency of the transformation. The protocol adopted by us in order to describe the physical and chemical processes described in this thesis seems to be reasonable. All the obtained results confirms that our methodology works good and the future perspectives are to apply the same protocol to another kind of molecular motor, artificial or natural, in order to describe it and obtain informations which cannot be obtainable experimentally and it also can be of great help during the design stage of them in order to obtain real world system with the wanted performances. In the second part of the work we have used the same approach of the first (problem described by PDEs, solutions and simulations with the FEM method, extraction and numerical elaboration of results) in order to describe one of the most studied and interesting phenomena that relies the life, the generation and propagation of nervous signals between excitable cells. This phenomena can be successfully described by the Hodgkin and Huxley mathematical model (Chapter 5) which can describe the generation and propagation of action potentials along the membrane of the active cells/fibres. We have thus applied all our know-how on this problem firstly validating it; in the first half of the Chapter 6 a huge work of validation is showed. We have simulate every kind of real-world experiment which can be done over an active cell such as the Voltage and Space clamp with results which matches perfectly the experimental evidences. Done this we gone over the monodimensionality of the HH model; coupling the solutions from it with the equation describing the space propagation of both electric fields and potentials allowed us to describe how an action potentials can a↵ect the space surrounding the active cell (second half of Chapter 6). With this results is finally


possible to exit the cell membrane and to go over the limits of the HH model and study a huge number of side-e↵ects related to the extracellular action potentials. The creation of successfully bidimensional models will lead to a huge quantity of future studies and researches. We cannot only simulate the extracellular action potential; one of the next steps can be to simulate and model the extracellular stimulation of active fibre (i.e. functional stimulation) and from this, model the ephaptic coupling of two (or more) fibres; one fibre acts a source for the extracellular potential which will act as a stimulator for the adjacent fibre(s). Is also possible to simulate a huge quantity of di↵erent extracellular media in order to understand which can be the best in order to improve such non-synaptic communications. Moreover our results and our model have found an immediate application inside an European project in which our research group in involved, the iOne-FP7 [2]. The target of this project is to develop and test for the first time an implantable device (Active Multifunctional Implantable Devices, AMIDs) which can lead to treatment of Spinal Cord Injury (SCI). In such a project our research group enters as a simulation support to the experimental tests and evidences and the correct computation of the extracellular action potentials, their range of action, their form and their intensity is of extreme importance for this project which can, in future, change the life of thousands peoples.

References [1] L. Moro, E. Bakalis, M. Di Giosia, M. Calvaresi, and F. Zerbetto, Operations and thermodynamics of an artificial rotary molecular motor, (2013), submitted. [2] http : //ione