Table of Contents Author Guidelines Submit a Manuscript
Advances in Mathematical Physics
Volume 2019, Article ID 7687643, 12 pages
Research Article

On Oscillations and Noise in Multicomponent Adsorption: The Nature of Multiple Stationary States

1Centre of Microelectronic Technologies, Institute of Chemistry, Technology and Metallurgy, University of Belgrade, Njegoševa 12, Belgrade, Serbia
2Faculty of Physical Chemistry, University of Belgrade, Studentski trg 12-16, 11000 Belgrade, Serbia

Correspondence should be addressed to Olga M. Jakšić;

Received 7 August 2018; Revised 25 November 2018; Accepted 16 December 2018; Published 1 January 2019

Academic Editor: Carlo Bianca

Copyright © 2019 Olga M. Jakšić et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Starting from the fact that monocomponent adsorption, whether modeled by Lagergren or nonlinear Riccati equation, does not sustain oscillations, we speculate about the nature of multiple steady state states in multicomponent adsorption with second-order kinetics and about the possibility that multicomponent adsorption might exhibit oscillating behavior, in order to provide a tool for better discerning possible oscillations from inevitable fluctuations in experimental results or a tool for a better control of adsorption process far from equilibrium. We perform an analysis of stability of binary adsorption with second-order kinetics in multiple ways. We address perturbations around the steady state analytically, first in a classical way, then by introducing Langevin forces and analyzing the reaction flux and cross-correlations, then by applying the stochastic chemical master equation approach, and finally, numerically, by using stochastic simulation algorithms. Our results show that stationary states in this model are stable nodes. Hence, experimental results with purported oscillations in response should be addressed from the point of view of fluctuations and noise analysis.

1. Introduction

Adsorption processes are surface phenomena that are well understood, explored, and employed in various scientific and industrial fields [14]. Some mathematical models used for the interpretation of adsorption kinetics are phenomenological, obtained through analyzing experimental data, while some are derived analytically. Either way, mathematically, adsorption can be modeled either as a process with pseudo-first-order kinetics (described by linear first-order differential equations, first proposed by Lagergren [5]) or as a process with second-order kinetics (described by nonlinear Riccati equations [6].) However, every mathematical model is a simplification of what is observed in real life experiments, and real life experiments are sometimes hard to interpret due to fluctuations and noise.

Natural systems are genuinely stochastic in nature. There are inherent intrinsic noise sources, extrinsic fluctuations originating from the surrounding equipment, and signal artifacts that also affect observed results. If deviations from ideal response in practical applications are reproducible, with some sort of purported regularity, we need a proper tool for interpretation of their origin with respect to difference between the oscillations which might emerge if proper conditions are met and the fluctuations which are inevitable. Oscillations and noise may relate in many different ways. If present in real life systems, they coexist. Sometimes small fluctuations and noise may invoke escalating oscillations in highly unstable systems, and sometimes oscillations are small and latently hidden in the noise, but sometimes noise may mimic oscillations and it is hard to deduce the origin of temporal deviations from the steady state, especially when the number of the samples is not very high. Stability is an important issue in many applications, as well as the appropriate modeling of the process under investigation.

Stability and chemical oscillations were analyzed in a classical way by Prigogine and coworkers in the early sixties. In [7] Nicolis and Portnow addressed the general criteria about the type of chemical reactions exhibiting sustained oscillations, available even before. Although these methods are still employed in practical situations [8], new approaches emerged with time. Very complex systems can be analyzed with greater ease numerically rather than analytically, and algorithms for simulations of systems, modeled as either deterministic or stochastic ones, enabled a more informative insight into specific reactions observed. Gillespie and colleagues proposed algorithms [916] that are frequently applied and also implemented into a standalone software solution for multiple platforms (StochKit2). These algorithms for simulation of stochastic reaction kinetics are based on the master equation for the probabilities that the system has a given composition as a function of time, i.e., chemical master equation (CME), the mandatory method in literature on stochastic analysis [17]. However, despite the completeness of the insight into the system one could gain by knowing any given composition at any instant, it is seldom the case that the CME can be solved analytically. Then, it is often the case that Langevin approach is applied, where a known deterministic model of the system is augmented with stochastic terms representing intrinsic noise sources and random fluctuations, with the main issue being the modeling of variances of the Gaussian zero mean added terms called the Langevin forces.

Apart from the focus on the reactants and products in reactions, there is also an approach where the reaction flux, the chemical reaction rate of the flow of reactants and products through the whole reaction network, is in the focus of the analysis, analogously to the analysis of generalized flux and its conjugated generalized force in the framework of nonlinear thermodynamics. The flux balance analysis approach enables the prediction of the growth of an organism in population dynamics or prediction of the rate of production of biotechnologically or chemically important metabolites and chemicals. The stochastic treatment of the flux is even more insightful because it reveals statistics about species creation and depletion [18]. Works of Bianca and Lemarchand on the reaction fluxes of the steady state and the oscillating biological and chemical systems showed how the relation between the reaction flux and the cross-correlation functions of the concentration fluctuations may be utilized in the analysis of biological and chemical systems [1921]. This approach has been so far the most appropriate for the analysis in this paper since it allows for addressing both oscillations and fluctuations at the same time.

Our motivation to contribute to this field emerged from the contradictory interpretations of the experimental results, where fluctuations can be misunderstood for oscillations and vice versa, and from the fact that adsorption is important beyond laboratories, and it is also an industrial process used in chemical, petrochemical, oil-refining, pharmaceutical, and related industries. In those complex systems, due to nonlinearities, small oscillations may escalate enough to endanger the whole operation of the chemical plant so additional investments in improvements of the automatized control of the system are necessary and economically justified [22]. Despite the fact that the phenomenon of oscillating adsorption process has been reported in praxis recently [23], we present here our doubts on the possibility of sustained oscillations in monolayer adsorption process itself, in spite of intrinsic nonlinearities in the model of the process [24].

We perform a theoretical analysis of the nature of multiple steady states in multicomponent monolayer adsorption with second-order kinetics in order to provide arguments for proper interpretation of seemingly repeatable deviations from expected values in experimental results. We perform the analysis in multiple ways. We addressed perturbations around the steady state analytically, first in a classical way, then by introducing Langevin forces and analyzing reaction flux and cross-correlations, then by the use of stochastic chemical master equation approach, and finally, numerically, by the use of stochastic simulation algorithms.

2. Materials and Methods

In order to analyze multicomponent adsorption using various approaches, we form the stoichiometric equation. In chemisorption, there is a chemical reaction between an adsorbed particle and an adsorption center, whereas it is not the case with physisorption, where no chemical bonds are made. However, in both cases adsorption is considered here as the process that is well described by the same mathematical formalism. The stoichiometric equation valid for general multicomponent adsorption of different adsorbate species in gas/liquid phase competing for free adsorption sites on the same surface, reversibly forming adsorbed particles , without dissociation or mutual interaction, iswhere and are adsorption and desorption rate constants, respectively. We assume a homogeneous surface with equivalent adsorption sites with no dissociative adsorption. Let us denote the total number of different species in the overall system by . In any given moment, of them are adsorbed whereas are in gaseous phase. There are no interactions between species in the gas phase. Let be the overall number of adsorption sites on the surface. There is a conservation relation that says that at any instant the sum of number of overall adsorbed molecules of all species on the surface and the number of free adsorption sites must equal . Considering all that, the set of deterministic differential equations that correspond to the stoichiometric equation (1) is as follows.These relations are also valid when the system reaches steady state. In the steady state the number of adsorbed molecules of each mixture component on the surface is constant and the adsorption rates equal the desorption rates. As a matter of fact, adsorption and desorption of molecules keep on occurring, making topological changes of adsorbate molecules on the surface but keeping the mean percentual surface coverage and mean number of adsorbed particles constant. These perturbations to a steady state will be addressed later. Let us now denote the numbers of adsorbate molecules on the surface for every component in the mixture in the steady state with . In the steady state the reaction flux of all relations in (2) equals zero. Equation (2) becomes a set of interconnected quadratic equations where each of the steady state numbers can be written in terms of the first.After replacing all relations from (3) into the first equation from (2), keeping the left hand side equal to zero, the steady state value appears to be a solution of a polynomial of the order of .Each solution to the polynomial , referring to the steady state number of adsorbate species denoted by one, generates the steady state numbers of all other adsorbate species, so it appears that the adsorption process of mixtures with adsorbate species may exhibit r+1 steady states in that particular system.

This result holds for an arbitrary mixture, even for the monocomponent adsorption, which is modeled by second-order Riccati differential equation. It is shown in [24] that only one of those two potential steady states obtained by mathematical calculus is practically possible; only one solution of the two is the steady state, with the other one being a mathematical artifact as it refers to a value greater than , the maximal number of adsorption centers on the surface, and also the maximal number of adsorbed molecules at any instant during the process of monolayer adsorption. The article [24] also addresses the criteria for modeling adsorption with approximate pseudo-first-order kinetic model while the article [28] addresses the stochastic analysis of adsorption of an arbitrary mixture for which the use of the approximate first-order kinetic modeling is appropriate. According to the classical criteria on stability [7], the linear closed systems are stable and do not exhibit sustained oscillations; the results from [28] are in agreement and also show that the mean value of the instantaneous number of adsorbed molecules of every species equals its deterministic value and that all the higher moments can be expressed in terms of the first moment (the mean value). Since the analysis of multicomponent adsorption stability for linear systems is resolved, the purpose of this work now is to investigate the nature of the steady states in adsorption of an arbitrary multicomponent mixture modeled as a process with the second-order kinetics only.

Inspired by the fact that any mixture may be approximated with fictive parameters referring to a single component gas, like using UDG option (user defined gas) or “one gas analyzer” option while operating binary gas analyzers [29, 30], we start the analysis of the possibility that the adsorption process of an arbitrary mixture may exhibit multiple steady states, with speculations referring to adsorption process of binary mixtures.

Stability is addressed firstly in a classical way, by analyzing the Jacobian of the Taylor expansion of the equation set (2), and then the procedure based on the Langevin approach is applied. Langevin forces and the relation between reaction flux and cross-correlation functions are modeled according to the procedures presented in [1921]. Then stochastic analysis, based on the chemical master equation, is applied. Final conclusion on the stability of adsorption of an arbitrary mixture is based on the fact that the analysis is general, with the only constraints being that the rate constants are positive real numbers and all numbers of molecules are integers. In the case kinetics and surface occupancy by any mixture can be modeled as if it were a single gas, just like when operating a binary gas analyzer, conclusions valid for the binary mixture can be generalized to more complex systems.

3. Results and Discussion

3.1. Classical Approach

The set of interdependent deterministic matrix Riccati differential equations valid for binary adsorption is as follows.Here, just like in general solution (4a) and (4b), the steady states, and , are interdependent. There are three possible sets of steady state numbers of adsorbate molecules in adsorption of a binary mixture. They can be calculated from the following equations.We seek the solutions in the form of positive real numbers smaller than or equal to the maximal possible number of the adsorbed particles, , keeping all other solutions neglected as a mathematical result without physical meaning. We start the analysis of the stability of steady states considering small perturbation of the system from a steady state. The kinetic equations set (5a) and (5b) written in the form of Taylor expansion around the steady state is as follows.The functions and equal zero in a steady state and the higher order terms may be omitted if perturbation from the steady state is small, so equation (8a) and (8b), written in the matrix form, is as follows. is the column matrix of the elements N1 and N2, and the expressions that correspond to the elements of the Jacobian are as follows.The solution to (9a) and (9b) implies an exponential regression to/from steady state or equilibrium.Here, are constants related to the initial perturbation from the steady state, while the time constants written in terms of the Jacobian are as follows.Bearing in mind that the Jacobian refers to the steady state, we analyze now the nature of the time constants. In order to determine the signs of , and of the discriminant in expression for lambda, we investigate the signs of , , , and first.

If we employ relations valid for steady states, we transform the expressions for and into the following.In monolayer adsorption the instantaneous number of adsorbed particles of any species can never be greater than the initial number of particles, and the instantaneous number of all occupied adsorption sites on the surface can never exceed , soand, hence, the Jacobian trace is always negative, and elements are always negative, and the discriminant in the expression for lambda is always positive.

Slightly more cumbersome algebra shows that the Jacobian determinant is always positive, which means that the discriminant root in (11) is always less than the Jacobian trace. Both of the time constants (the eigenvalues of the Jacobian matrix) are always real negative numbers. That means that the steady state is a stable node and oscillations cannot be sustained. All departures from the only practically feasible steady state end up with exponential regression back towards that one and the only steady state without the chance to reach the state of a Hopf bifurcation point under any circumstances.

By analyzing the stability of binary adsorption with this procedure, we came to conclusion about the nature of the steady states without calculating or knowing their exact value. We analyzed the deviations from the steady state and the relations valid for the steady state, but the value of the steady state was not needed for the application of criteria in this procedure. This is important because the analysis of the cubic polynomial is not needed.

So far we did not address the stochastic nature of the system, its intrinsic fluctuations, and causes or frequencies of the deviations from the steady state.

3.2. Langevin Forces and Correlation Functions

In this approach we focus on the part of the system we are interested in as if it were purely deterministic, and we consider the rest as a statistical bath [17]. We do not consider the Taylor expansion of functions in kinetic equations now; instead we rewrite the original kinetic equations replacing the instantaneous number of adsorbed molecules with the sum of its steady state value and the instantaneous deviation from the steady state . As before, the indices refer to a particular species in the binary mixture.We assume that the departures from the steady state are small so it is convenient to neglect all the products resulting from their multiplication. We also neglect all terms related to the reaction flux, which does not comprise infinitesimally small variables but vanishes at the steady state. However, we add the Langevin forces, random terms to interpret randomness in departures and as a compensation for all of the neglected terms. The added terms are zero mean Gaussian noise terms and whose variances will be addressed later in this section. By doing so we linearized the equation set (15a) and (15b). The constants in (16a) and (16b) are defined by the following expressions.In order to investigate the stability of the steady state, we rewrite the equation (16a) and (16b) in the matrix form and compute the eigenvalues of the matrix .Now eigenvalues are as follows.Bearing in mind that expressions in (18) refer to the sorption rate constants, the overall numbers of molecules of reactant species N01 and N02, and the overall number of adsorption centers on the surface that are all real positive numbers, we see that the trace of the matrix is always negative as it never reaches zero, which means that a Hopf bifurcation can never occur. That is in accord with the stability analysis performed without Langevin linearization.

Langevin forces have zero mean but nonzero variance. In order to determine the variances of the Langevin forces, we focus on the stoichiometric equations and the reaction flux in equilibrium. In case the reactions are reversible, as they are in the model of adsorption, forward and backward fluxes are treated separately. In the case of binary adsorption we have the following.Although movements of molecules of different species are independent in the same manner as they are in any Brownian motion with Maxwell-Boltzmann distribution, on the surface they fight for the same area so we have the additional conservation equation that says that at any moment the sum of the number of adsorbed molecules and and the number of free adsorption sites must equal the overall number of adsorption sites on the surface, . Additionally, for each of both species the initial number of molecules in the system, N0, must equal the sum of the number of its adsorbed molecules and the number of the molecules in the gas phase .

For every stoichiometric equation and every variable we are interested in, we define one stochastic term; let us denote it by multiplied by the square root of adsorption or desorption rate in the steady state. Then we form the Langevin force for every stochastic variable, taking into account all stoichiometric equations in which it appears.

Langevin forces in (16a) and (16b) are given by the following.All stochastic terms are zero mean Gaussian white noise terms with unit variance, mutually uncorrelated, with their cross-correlations being zero. Hence, the variances of the Langevin forces arewhereIn order to gain a better insight into the adsorption process and how it behaves at very short times and high frequencies, we focus now on the correlation functions. We use the procedure in [1921] and determine the correlation functions after the change of the basis by the use of the following matrix. The matrix relates and with new coordinates.After finding the correlation functions of the new coordinates and using the change of the basis, in the long time limit when all transients are gone and only the dependence on the time difference remains, then the difference of cross-correlations of the deviations to the steady state is given by the following.After multiple rearrangements, by using all expressions valid for the matrix elements, the variances of the Langevin forces, and time constants in (26) and by using all relations valid for the steady state, the expression for the difference of cross-correlation functions is obtained as follows.The reaction flux, related to the original equation set (5a) and (5b), before linearization, vanishes at equilibrium and equalsAlthough the Langevin approach is widely used in practical applications, it does not give the full insight into all statistical parameters of the process as a fully stochastic approach by the use of the chemical master equation would do.

3.3. Stochastic Approach: Chemical Master Equation

Instead of investigating deterministic equations in the vicinity of a steady state or stochastic Langevin differential equations with added random white noise, here we analyze adsorption as a stochastic process and work with probabilities that instantaneous numbers of adsorbate molecules on the surface may have a specified value. By doing so, we apply the conventional procedure of stochastic analysis and assume that the time span for our system observations, , is short enough, so that only one transition between the neighboring states may occur. That means that the system can reach the state with N1 adsorbed molecules of the first species and N2 adsorbed molecules of the second species in five ways: one molecule of the species denoted by 1 has been adsorbed, one molecule of the species denoted by 2 has been adsorbed, one molecule of the species denoted by 1 has been desorbed, one molecule of the species denoted by 2 has been desorbed, and no transition occurred. The probability that the system will be in a state with N1 adsorbed molecules of the first gas species and N2 adsorbed molecules of the second gas species in the moment is the sum of these five probabilities.The probability that no transition occurs is as follows.The transition probabilities among neighboring states are proportional to the instantaneous sorption rates, soIn the short time limit, after rearranging terms in (29), we obtain the differential equation for the probabilities, i.e., the chemical master equation.In order to find the solution to the chemical master equation governing the adsorption process of binary gas mixture, we need a bivariant probability generating function. Its first derivative with respect to , calculated for s1 = 1 and s2 = 1, is then used for the calculations of the first moment of the number of adsorbed molecules for the first gas species.Its first derivative with respect to s2, calculated for s1 = 1 and s2 = 1, equals the first moment of the number of adsorbed molecules for the second gas species.Its second derivative with respect to s1 and s2 is as follows.Sums go to infinity but their terms are zero beyond physical limits. In order to find the solution to CME (32) we first multiply it by and then by and then rearrange its members grouping together the terms that form the probability generating function (PGF) or the derivatives of PGF. By doing so, we transform CME, differential equation on probabilities (32), into a partial differential equation on PGF, F(s1,s2,t). Once solved, the partial differential equation on F(s1,s2,t), would reveal all statistical data on the process with all higher moments.

Since it is a cumbersome and possibly unfeasible task, with no great importance for the applications where knowing of the first moments only is crucial, it is often much faster to seek for the solution by simulations or the analysis of the moments only. Here, we will address it using simulation algorithms based on the chemical master equation.

3.4. Stochastic Approach: Simulation Algorithms

Our routines for visualization of stochastic transitions in multicomponent adsorption of molecules are based on stochastic simulation algorithms (SSA) developed by Gillespie and colleagues [916]. The routines are part of our custom, in-house designed, permanently growing software solution. Features of its current version, ADmoND, written for MathWorks MATLAB release 2013, will be briefly outlined in next section and the whole package is available at Mendeley Data repository [31].

The propensity functions, used for modeling the probabilities of transitions between states and for the modeling of the time sequence of random moments for every transition, are the same as in the chemical master equation (32). Here, the mean time between transitions of the molecules of each species in the mixture is a random variable with the exponential distribution whose parameters are related to propensity function, i.e., sorption rate of the relevant transition (adsorption or desorption of the first or second component in the mixture). The exponential distribution is designed using in-built function for uniform distribution in MATLAB and inverse CDF sampling technique method [32]. After generating the sequences of time and the instantaneous numbers of adsorbed molecules for each component in the mixture, it is easy to obtain the time evolutions of the number of adsorbed molecules for every component in the mixture, the instantaneous changes in the phase space, and the correlation functions.

Here we give results obtained by the simulation algorithm valid for binary adsorption of the mixture of Sarin and Sulphur Mustard on graphene surface. Table 1 lists parameters of the mixture. The sorption rate constants are calculated by the procedure from [24], and the desorption energy and the surface density of adsorbed particles, needed for the calculations of the sorption rates, are taken from the online Harvard Dataverse repository [25] and from [26, 27]. Assumed volume of the reaction chamber is 1 dm3, temperature is 282 K, and the number of adsorption centers on the surface, , is 38.17.1014 (the surface density of graphene atoms being 38.17.1018 and surface area being 1 cm2). Figure 1(a) shows phase space assuming that, at the initial moment, the surface was clean with no adsorbed amount. Figure 1(b) shows the corresponding cross-correlation function. Due to very high desorption rate constants (Table 1) and high number of adsorption sites on the surface, 1 million transitions were not enough to show the reaching of equilibrium. The corresponding correlation function is symmetric.

Table 1: Parameters used in simulations of binary adsorption on graphene in Figure 1.
Figure 1: Phase space (a) and the corresponding cross-correlation function (b) of numbers of adsorbed molecules of Sarin and Sulphur Mustard on graphene assuming clean surface in the beginning.

Figure 2 shows simulations of adsorption of fictive binary mixture assuming the number of adsorption centers on the surface, , is 38.17.105, adsorption rate constant and desorption rate constant of the first component 2.89 and 1.7, respectively, and adsorption rate constant and desorption rate constant of the second component 1.17 and 1.57, respectively, assuming the initial condition is the clean surface with no adsorbed molecules. Other parameters are the same as for Figure 1. On the left, up, and down, on time evolutions of numbers of adsorbed molecules, we see either that reaching the equilibrium was fast and that there was just one transient or that we need more than 10 million transitions in simulations in order to see vanishing of the other transient. In order to see exact equilibrium levels, we would have to check with the analytical solutions. On the right, we see that every simulation led to the same node and that the cross-correlation functions were symmetrical.

Figure 2: Time evolutions of the number of adsorbed molecules (left), phase space (up right), and the corresponding cross-correlation function (down right), fictive mixture, clean surface at the beginning.

Figure 3 shows results for the same fictive mixture as Figure 2, with the only difference being the starting moment of observation, which is now random. Phase space, shown on the upper diagram on the right side, shows that all different starting points are leading to the same node. The diagram on the lower right part in Figure 3 shows that the cross-correlation function is not ideally symmetric in this case. Time evolutions in Figure 3 (top and bottom diagram on the left) demonstrate now two changes in the process kinetics. Due to our analysis we now know that there are no two consecutive steady states, and change in kinetics is the consequence of the vanishing of multiple transients (for binary mixture the exact solution would give at least two transients). Exact solution to multicomponent adsorption modeled with linear Lagergren model, published in article [28], states that in linear systems the number of transients equals the number of components in the mixture. The exact solution to nonlinear systems would reveal that there are more transients than there are components in the mixture but they may not be always visible in the observed time and amplitude span. The quality of results also depends on the quality and features of used software.

Figure 3: Time evolutions of the number of adsorbed molecules (left), phase space (up right), and the corresponding cross-correlation function (down right), fictive mixture, random initial moment.
3.5. Algorithms and ADmoND Software for Adsorption Analysis

The software package ADmoND we used is aimed at the analysis of monolayer adsorption and desorption processes in novel devices. It has been developed for the software platform MathWorks MATLAB, 2013a release. ADmoND can be used as a MATLAB toolbox or MATLAB application with graphical user interface (GUI) that is adapted MATLAB GUI for figure handling, extended with additional tabs (System, Evolutions, and Equilibrium). All calculations are based upon analytically derived expressions, implemented in mathematical forms optimized to avoid or minimize numerical error propagation. Theoretical background is based on well-known laws or new original solutions published in peer-reviewed journals. The full software package is openly available and can be downloaded from Mendeley Data repository [31].

The 'System' tab is aimed at handling parameters saved outside MATLAB and interactive setting of system parameters (temperature, adsorbent area, reaction chamber volume, and composition of an adsorbate mixture). Figure 4, generated with ADmoND, illustrates adsorption of binary gas mixture on a solid surface.

Figure 4: Competitive adsorption between particles of two species for equivalent adsorption sites on homogeneous surface.

The 'Evolutions' tab of ADmoND is dedicated to calculations related to adsorption kinetics in time domain. Figure 5 illustrates possible time evolutions of monocomponent adsorption (green lines) obtained by stochastic simulation algorithm (SSA, inspired by [9], adapted for ADmoND case studies), their mean (magenta line), and the ideal deterministic solution (red dashed line).

Figure 5: Time evolution of monocomponent adsorption obtained by stochastic simulation algorithm (SSA, green lines), mean of SSA ensemble (magenta line), and deterministic solution (red dashed line).

The 'Equilibrium' tab is dedicated to calculations related to adsorption dynamics in equilibrium. As we saw, for every mixture of components, there are r+1 possible sets of numerical values, r+1 possible steady states. Each possible steady state is a vector of r numbers, where each number corresponds to the number of adsorbed molecules for every component in the mixture (). ADmoND finds r+1 roots of the polynomial of the r+1 order (candidates for the steady state of component 1) then computes all other possible steady states for all other components in the mixture. Afterwards, it performs inspection if exactly one set of solutions for the number of adsorbed particles of mixture components in steady state Ns is realistic and satisfies three conditions. The first condition is that the steady state values for the number of adsorbed molecules of particular components must be positive numbers, lower than NM. The second condition is that the sum of all steady state values for the number of adsorbed molecules of particular components on the surface may never exceed . The third condition is that the steady state values for the number of adsorbed molecules of each component in the mixture must be positive number, lower than the overall number of its molecules in the system, .

The theoretical proof that there is just one steady state, the stable node, ensures that there is just one set of solutions and thus enables the optimization of the code (the first set that is found is the only set; there is no need for inspection of all other sets).

3.6. Discussion of Results

The stability analysis, described above, helps in a better understanding of adsorption dynamics in equilibrium and a more accurate interpretation of results. In research span of almost a decade, in simulations of multicomponent adsorption, we did not encounter a case of practically feasible multiple steady states obtained by MATLAB; only one set of solutions satisfied required conditions (steady state adsorbate molecules on the surface are positive numbers with overall sum that is smaller than or equal to the maximal possible number of adsorbed particles, ). After thorough investigation, some complex solutions that appeared once proved to be the consequence of the numerical error propagation, which resulted in our additional efforts towards code improvements.

All of the methods we applied here for stability analysis in this paper result with the same conclusion: There is no combination of system parameters that ensures the appearance of sustained oscillations. On the contrary, the steady state is the stable node. Thus, if there is a seemingly oscillatory response in the experiment, the cause of oscillations is not the adsorption process (modeled with the nonlinear set of Riccati equations). Our conclusion is that the possible explanations are as follows:(i)There are no oscillations; noise is mimicking them.(ii)The reactions are more complex and are not well modeled by the set of Riccati equations.(iii)The oscillations emerged due to some other cause in the control loop of complex systems used in chemical, petrochemical, oil-refining, pharmaceutical, and related plants.

An example that supports this reasoning is presented in the paper [23], where self-oscillating adsorption–desorption of silver nanoparticles on the gel formed by incorporating graphene oxide into a poly (methacrylic acid)–polyethylene glycol copolymer film is reported. Here, adsorbent is a pH-responsive hydrogel which acts as a pH oscillator in a closed reaction and the oscillations in reversible adsorption–desorption of Ag nanoparticles on the patterned hydrogel surface are just response to pH oscillations.

4. Conclusions

The presented results prove that a binary adsorption process, moved from the equilibrium, goes through different transition states while reaching a new equilibrium, without experiencing oscillations. It is valid for all processes governed by the analyzed equations or their mathematical equivalents, i.e., for adsorption of gases but also for biological antibody-antigen kinetics or adsorption of macromolecules observed in commercial surface plasmon resonance-based instruments. If we extend this reasoning to mixtures with more than two components, we see that this result is also valid. It is valid whenever we can divide a mixture into a target gas and a fictive gas with user defined parameters (analogously to “user defined gas” option in commercial binary gas analyzers [30]), if the before mentioned assumptions are valid (monolayer adsorption, noninteracting species, and positive real sorption rates).

This result can be used for better understanding of seemingly oscillating experimental values, for better judgment of experimental setup, or for better control of adsorption systems far from equilibrium. Further efforts in research and development of adsorption based devices useful for operation of commercial surface plasmon resonance instruments or chemical plants in petrochemical, oil-refining, paper-producing, pharmaceutical, and related industries might be directed to the field of analysis of noisy response.

Data Availability

Data on desorption energy and surface density of adsorbed molecules of chemical warfare agents, needed for calculations of the sorption rate constants, used here for the demonstration of theoretical results, are available on the Harvard open source research repository at [25]. ADmoND software package written for the MathWorks MATLAB environment is available on open research repository Mendeley Data at [31].

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This work has been supported by the Serbian Ministry of Education and Science through projects TR 32008 and ON 172015. The authors acknowledge Dr. Annie Lemarchand, CNRS, Paris, France, for insightful correspondence regarding modeling the Langevin forces.


  1. S. Motoyuki, Adsorption Engineering, Elsevier, 1990.
  2. D. D. Do, Adsorption Analysis: Equilibria and Kinetics, vol. 2, Imperial College Press, 1998.
  3. W. J. Thomas and B. Crittenden, Adsorption Technology and Design, Elsevier Science & Technology Books, 1998.
  4. A. Da̧browski, “Adsorption—from theory to practice,” Advances in Colloid and Interface Science, vol. 93, no. 1–3, pp. 135–224, 2001. View at Publisher · View at Google Scholar · View at Scopus
  5. Y. Liu and L. Shen, “From Langmuir kinetics to first- and second-order rate equations for adsorption,” Langmuir, vol. 24, no. 20, pp. 11625–11630, 2008. View at Publisher · View at Google Scholar · View at Scopus
  6. J. Toth, Adsorption: Theory, Modeling, and Analysis, Surfactant, Marcel Dekker Incorporated, Surfactant, 2002.
  7. G. Nicolis and J. Portnow, “Chemical Oscillations,” Chemical Reviews, vol. 73, no. 4, pp. 365–384, 1973. View at Publisher · View at Google Scholar · View at Scopus
  8. V. Cupic, A. Ivanovic, and L. Kolar-Anic, “Modeling of the Complex Nonlinear Processes: Detrmination of the Instability Region by the Stoichiometric Network Analysis,” in Mathematical Modeling, C. Brennan, Ed., pp. 111–178, Nova Science Publishers, 2011. View at Google Scholar
  9. D. T. Gillespie, A. Hellander, and L. R. Petzold, “Perspective: Stochastic algorithms for chemical kinetics,” The Journal of Chemical Physics, vol. 138, no. 17, Article ID 170901, 2013. View at Publisher · View at Google Scholar · View at Scopus
  10. D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” The Journal of Physical Chemistry C, vol. 81, no. 25, pp. 2340–2361, 1977. View at Publisher · View at Google Scholar · View at Scopus
  11. D. T. Gillepsie, “A general method for numerically simulating the stochastic time evolution of coupled chemical reactions,” Journal of Computational Physics, vol. 22, no. 4, pp. 403–434, 1976. View at Publisher · View at Google Scholar · View at Scopus
  12. D. T. Gillespie, “Approximate accelerated stochastic simulation of chemically reacting systems,” The Journal of Chemical Physics, vol. 115, no. 4, pp. 1716–1733, 2001. View at Publisher · View at Google Scholar · View at Scopus
  13. S. Lampoudi, D. T. Gillespie, and L. R. Petzold, “The multinomial simulation algorithm for discrete stochastic simulation of reaction-diffusion systems,” The Journal of Chemical Physics, vol. 130, no. 9, Article ID 094104, 2009. View at Publisher · View at Google Scholar · View at Scopus
  14. D. T. Gillespie, “The Chemical Langevin equation,” The Journal of Chemical Physics, vol. 113, no. 1, pp. 297–306, 2000. View at Publisher · View at Google Scholar · View at Scopus
  15. D. T. Gillespie, “The chemical Langevin and Fokker-Planck equations for the reversible isomerization reaction,” The Journal of Physical Chemistry A, vol. 106, no. 20, pp. 5063–5071, 2002. View at Publisher · View at Google Scholar · View at Scopus
  16. D. T. Gillespie, “Stochastic Simulation of Chemical Kinetics,” Annual Review of Physical Chemistry, vol. 58, no. 1, pp. 35–55, 2007. View at Publisher · View at Google Scholar
  17. C. W. Gardiner, Handbook of Stochastic Methods: For Physics Chemistry and the Natural Sciences, Springer, New York, NY, USA, 1985. View at MathSciNet
  18. O. Kahramanoǧullari and J. F. Lynch, “Stochastic flux analysis of chemical reaction networks,” BMC Systems Biology, vol. 7, article no. 133, 2013. View at Publisher · View at Google Scholar · View at Scopus
  19. C. Bianca and A. Lemarchand, “Evaluation of reaction fluxes in stationary and oscillating far-from-equilibrium biological systems,” Physica A: Statistical Mechanics and its Applications, vol. 438, pp. 1–16, 2015. View at Publisher · View at Google Scholar · View at MathSciNet
  20. C. Bianca and A. Lemarchand, “Determination of reaction flux from concentration fluctuations near a Hopf bifurcation,” The Journal of Chemical Physics, vol. 141, no. 14, Article ID 144102, 2014. View at Publisher · View at Google Scholar · View at Scopus
  21. C. Bianca and A. Lemarchand, “Temporal cross-correlation asymmetry and departure from equilibrium in a bistable chemical system,” The Journal of Chemical Physics, vol. 140, no. 22, 2014. View at Google Scholar · View at Scopus
  22. Online Oscillation Detection, PiControl Solutions Company, 2015,
  23. J. H. Jang, M. Orbán, S. Wang, and D. S. Huh, “Adsorption-desorption oscillations of nanoparticles on a honeycomb-patterned pH-responsive hydrogel surface in a closed reaction system,” Physical Chemistry Chemical Physics, vol. 16, no. 46, pp. 25296–25305, 2014. View at Publisher · View at Google Scholar · View at Scopus
  24. O. M. Jakšić, Ž. D. Čupić, Z. S. Jakšić, D. V. Randjelović, and L. Z. Kolar-Anić, “Monolayer gas adsorption in plasmonic sensors: Comparative analysis of kinetic models,” Russian Journal of Physical Chemistry A, vol. 87, no. 13, pp. 2134–2139, 2013. View at Publisher · View at Google Scholar · View at Scopus
  25. O. Jakšić, Parameters for adsorption of (CWA) chemical warfare agents and CWA simulants, Harvard Dataverse, 2018.
  26. B. N. Papas, I. D. Petsalakis, G. Theodorakopoulos, and J. L. Whitten, “CI and DFT studies of the adsorption of the nerve agent sarin on surfaces,” The Journal of Physical Chemistry C, vol. 118, no. 40, pp. 23042–23048, 2014. View at Publisher · View at Google Scholar · View at Scopus
  27. L. Ebrahimi, A. Khanlarkhani, M. R. Vaezi, and M. Babri, “Molecular dynamics simulation of the interface between sulfur mustard and graphene,” Computational Materials Science, vol. 152, pp. 355–360, 2018. View at Publisher · View at Google Scholar
  28. O. M. Jakšić, Z. S. Jakšić, Ž. D. Čupić, D. V. Randjelović, and L. Z. Kolar-Anić, “Fluctuations in transient response of adsorption-based plasmonic sensors,” Sensors and Actuators B: Chemical, vol. 190, pp. 419–428, 2014. View at Publisher · View at Google Scholar · View at Scopus
  29. User Manual BGA244 Binary Gas Analyzer, Stanford Research System,
  30. Report, “Binary Gas Analyzer T750 Portable Calibrator, Teledyne Api, 2015,
  31. O. Jakšić, “ADmoND: MathWorks Matlab Package for simulation of monolayer adsorption processes in nano devices,” Mendeley Data, vol. 1, 2018. View at Google Scholar
  32. S. M. Ross, Simulation, Elsevier Science, 2006.