Mathematical Modeling in Secondary Chemistry Education: Chromatography

The rapid advance in information technology requires further developments in all areas of education. In this context, one should think about going beyond the use of digital media for the mere presentation of scientific content. Interactive computer simulations allow quasi-experimental investigations of scientific phenomena but for students they usually remain black-box approaches. For a deeper understanding of phenomena, it is desirable to go one step further and set up computer codes based on a given microscopic model as part of the chemical education. Such approach allows teaching the scientific topic in more depth, fosters the awareness of the relevance of mathematics and computing in chemistry, and lastly supports the self-directed investigation of a scientific phenomenon. In addition, it gives students the opportunity to learn in general about modelling which has become an important contribution to chemistry and other natural and engineering sciences. Here we discuss basic chromatography with a simplistic stochastic simulation method suitable for upper secondary education. In addition, the analytical solution of the processes is given at the level of secondary mathematics. Chromatography itself is potentially treated in secondary education at various levels from paper chromatography to gas chromatography. This general knowledge makes it more accessible to students as a subject for deepening by modeling and simulation.


Introduction
The focus of this contribution is to provide a topic in secondary education that is suitable for introducing modelling and simulation in chemistry and is related to mathematics at the same level. The aims of this task are manifold. Besides introducing modelling and simulation it is focused on the deepening of chemical knowledge in chromatography and further related chemical topics.

Modelling and Simulation
Without doubt, computational methods are important in chemical sciences research. However, in the preparation of students in secondary education, modeling and simulation are less visible. There is a general proposal to teach what is called computational thinking in computer science. This idea has been proposed and developed several times in recent decades [1]. The most recent call for implementing computational thinking in K-12 education (kindergarten through grade 12) was by Wing [2]. Wing claims that computational thinking is of importance for everyone not just computer scientists. No matter whether students who likely will not enter computer science benefit from this way of thinking or not, it is important to give students in upper secondary education the opportunity to learn about the existence of such methods in general and in the context of chemistry. They should learn that experiments and theoretical methods complement each other and in this way get a more realistic view on chemistry.
Even though the focus here is rather on chemistry than on computational thinking itself, we can profit from those ideas. It is stated [1] that the most appropriate definition of computational thinking is given by Aho [3]: "thought processes involved in formulating problems so their solution can be represented as computational steps and algorithms". More generally, computational thinking can be regarded as "device independent computational process" [4] being the intersection set of electronic computing and biologic computing by humans. So computational thinking should rather be separated from developing computer codes only but may be regarded as procedural thinking [5]. We can use this approach to set up a computer code for a chemical problem by starting with a general procedure. Here we realize this by a board game with a game instruction acting as a procedure. This procedure can be used to set up a pseudo code, which is not yet a runnable computer code but a device independent computational process. Finally, we transfer the pseudo code to a real computer code that performs the simulation.

Chromatography
Chromatography can be introduced experimentally in secondary education at several levels. We can use paper chromatography already in lower secondary education, for example, for the decomposition of plant colors or felt pen colors. In upper secondary education, we can use gas chromatography with a simple column and a thermal or electric conductivity detector for decomposing lighter gas. Besides these practical applications, chromatography allows to treat several general phenomena in chemistry such as diffusion and molecular interactions.
Here we propose a stochastic simulation of chromatography feasible in high school. The resulting simulation code gives the students the opportunity to investigate the system with their simulation program by variation of the input parameters. Such modeling allows relating the particle model of matter to measurable macroscopic properties by statistical methods. While the behavior of each single particle is randomized and hence not predictable, the resulting macroscopic behavior follows the laws of stochastics. In this way, students learn from a simple example how a simulation works in principle beyond their perception of a black box. In case of stochastic simulations, it is additionally possible to introduce the relation between the molecular interaction energy and stochastics. Of course, there are many other approaches such as molecular dynamics simulations, mesoscale simulations, or continuum simulations, which should be mentioned to students at least to avoid the impression that all simulations are based on random numbers.
The fundamental molecular theory of chromatography goes back to Giddings and Eyring in 1955 [6]. Felinger has published a recent review in this field [7]. As a process being related to diffusion and interactions, several authors have modeled chromatography by stochastic simulations [8,9,10,11]. This includes also investigations of the underlying diffusion with the Ehrenfest urn model [12] rather on a fundamental level [13,14].
In educational science, chromatography is treated also for a long time. An early introduction goes back to 1969 [15]. Over the years several games have been proposed to introduce chromatography at the upper secondary and introductory level [16,17,18,19,20,21,22].
The educational approach here is to employ chromatography as an example for a theoretical investigation. Similarto scientific research approaches the process is divided into several detailed contributions. In order to relate microscopic processes to macroscopic phenomena one usually begins with the simplest possible model for a certain effect. Once this model is understood, one may successively add additional microscopic effects and observes the resulting macroscopic behavior of the system. This allows relating changes in macroscopic observations to the microscopic contributions. Here, within chromatography there are five steps: •Diffusion •Chromatography without stationary phase •Chromatography with stationary phase •Reconstruction of a chromatogram on the time axis •Deviation from the linear adsorption isotherm As we will see, each of these steps leads to a characteristic change in the chromatogram. Still the approach used here is relatively simple and covers only the basics. It is not meant to model experimental data quantitatively. A quantitative predictive simulation of chromatography requires to model molecular interactions in detail and much more advanced simulation techniques as, for example, described in a work on reversed-phase liquid chromatography [23].
Depending on the local curricula, the simulation of chromatography in chemistry also allows crosslinking to mathematics. This focuses on the binomial and the Gaussian distribution, which are obtained in the simulations here.

Methods
The process of setting up a simulation for chromatography requires several steps. Here we begin with the diffusion of a substance from a given starting position with fixed center of gravity of the dispersed substance. This means that the substance is not moving with a mobile carrier substance. The second step involves moving the center of gravity along a (virtual) tube. The third step is the inclusion of a stationary phase interacting with the substance. This leads to different retention times and allows separating two substances. Then, an experimental chromatogram is reproduced on the time axis. Finally, peak tailing is considered. Besides the scientific systematics this stepwise introduction allows deciding how detailed the topic is treated in a specific course which may depend on the educational level and capability of the students.

Diffusion
Diffusion was first described as Brownian motion [24] which is the random movement of macroscopic particles in liquids and gases. Fick set up two equations for modelling diffusion [25]. The first law of Fick describes a steady state diffusion with a constant concentration gradient while the second law accounts for changes in the concentration gradient. Einstein obtained the Gaussian distribution as solution for the Brownian motion. In this context, he also obtained the relation of the diffusion coefficient to the mean square displacement of the molecules [26].
Here, we start with the simulation of a one-dimensional diffusion of a substance around its center of gravity at position = 0 known as random walk model. To introduce the simulation method and to get a simple starting point for setting up a computer algorithm we begin with a board game. This game has been successfully employed here in a math course as an alternative to the Galton board [27] for the introduction of the binomial distribution. The only mathematical difference is the transformation from the play stone position to the index of the binomial distribution as discussed below.
In Figure 1a the initial setup for 12 stones representing 12 molecules diffusing along the -axis is shown. The stones are either numbered or in numbered lanes. All stones start at position = 0 . The -coordinate is numbered at the bottom of the board. Now we chose at random one of the 12 stones with equal probability of 1/12. This can be realized by a dodecaeder dice or by a random number generator app for mobile phones 1 . An app also allows using more than just 12 stones which it advisable in order to get a more meaningful distribution function. Once a stone is chosen we have to decide whether to move it into positive or negative direction. For this, a coin may be used. Then we move the chosen stone in the chosen direction and repeat the whole process many times. It is important to find a good end point because over time stones will leave the game board. For example, it turned out to be a good choice for 12 stones to stop when the first stone reaches = −4 or = 4 . Figure 1b shows a possible outcome with the distribution at the top. Due to the small number of stones, the distribution function is statistically not very meaningful but we can average the data of all student groups in class to obtain a more realistic distribution. In this context is important that all averaged data are for the same number of steps. One scope of game board besides the activity is to introduce the procedure for setting up a code as discussed below.  With a worksheet for a stochastic tree as depicted in Figure 2, the students can approach the distribution function by counting the number of paths that end at a given position. For example, after four time steps ( = 4) there is only one way to reach = −4 namely moving four times to the left. To reach = −2 there are four paths, depending at which step the move to the right is done. The number of paths to the different positions turn out to be 1 4 6 4 1 which are apparently the binomial coefficients that may be familiar to students from the Pascal triangle as well. To verify this we can do the same analysis for = 3 and = 2 . The probability for each single path is the product of the probabilities along the path. Here for the single step probability of = 1/2 and four steps it follows a probability 1/16 for each single path. Multiplying the number of paths with the single path probability 0,0625 gives the binomial distribution depicted in Figure 3. The general equation for the binomial distribution is given by: Here we have to apply the transformation from the positions of the stones along the -axis to the index of the binomial distribution with = ( + )/2 (see Figure 3).
In the next step we link the mathematical model to the physical-chemical properties. According to Fick's laws [25] the diffusion coefficient represents the magnitude of diffusion. It is related to the mean square displacement ∆ 2 by [26]: This specific equation is valid for one-dimensional diffusion as considered here. The mean square displacement can be identified by the standard deviation of the binomial distribution: Generally, in equation 2 and 3 the variable is the time in seconds. Here we may identify the real time with the number of steps , which justifies the use of the same symbol. Doing so and inserting equation 4 yields here = 1/8. A rigorous derivation of the standard deviation for the one-dimensional random walk model is given in the appendix. Since this is beyond secondary mathematics one may make the solution plausible to students based on their knowledge of the binomial distribution and its expression for the expectation value and the standard deviation. For this approach we start with the equation for the regular standard deviation of the binomial distribution on the axis of its index : With the equation for the expectation value of the binomial distribution on the -axis = • , and with = 1/2 we obtain: Now we apply the transformation from the index of the binomial distribution to the -axis ( = 2 − ) to the expectation value. We may argue to neglect the term − because the standard deviation -being the width of the distribution -does not change when we shift the distribution function along the abscissa. Finally, we obtain: The expectation value on the scale of the binomial index is the average number of moves per stone = run / A . Here run is the number of moves and A is the number of molecules of type A. Hence, we may compare the accuracy of the game results to the analytical solution: For this, the standard deviation of the board game data can be calculated using = 0 and the number of stones at position by: From students who are familiar with the binomial distribution from mathematics as for example in German upper secondary education, one can expect the knowledge of equations 1, 4, 8, and the equation for the expectation value. The mean square displacement has to be introduced at this level. The transformation between -and -axes can be elaborated by students if they know equation 1 and have evaluated the tree in Figure 2.
However, the accuracy of the game is not sufficient. It requires much more stones and steps to obtain a reasonable distribution. Therefore, we set up a computer algorithm that simulates the game and allows enlarging the system significantly. Due to the simple game rules, it is relatively straightforward to set up the corresponding computer algorithm. Such a code based on the game instruction is listed in Table 1(Supporting Information). Furthermore, a runnable code written for an educational Python application is provided in the Supporting Information. With respect to computer science, the code contains several basic programming structures such as loops and conditional statements. During the calculation, the data are saved in a file for plotting and further analysis. In Figure 4 the results for a simulation with 10 4 stones for 2 • 10 6 and 8 • 10 6 moves are shown. The solid curves are calculations with the Gaussian distribution using the input parameter = � run / A . This is possible because the binomial distribution can be approximated by the Gaussian distribution according to De Moivre and Laplace for • • (1 − ) > 9. For = 0 it is given by:

Taylor Dispersion like-Chromatography
A difference between diffusion and chromatography is the transport of the substance with a mobile phase. Hence, the center of gravity is not fixed = 0 but moves with a certain velocity while the diffusion takes place around the center of gravity. If there is no stationary phase that interacts with the moving substance, we observe the socalled Taylor dispersion [28,29]. This method is used for example to measure diffusion coefficients of solutes in supercritical fluids [30]. For this, the solute is injected in a flow of supercritical carbon dioxide in an empty tube. In that application the diffusion coefficient is obtained from the widening of the peak. Due to the simplistic approach here, we omit any characteristics of the flow such as a laminar velocity profile, which are important within the Taylor dispersion method [28]. So we just have a process with similarity to the full Taylor dispersion. The corresponding board game as well as the computer algorithm for the Taylor dispersion-like simulation is even simpler than for local diffusion because the stones need to be moved only in one direction. The initial state is depicted in Figure 5. With dice or a random number generator app one of the stones is chosen with equal probability. This stone in then moved one step towards increasing x-value. The corresponding computer algorithm is listed in Table 2(Supporting Information). Here, the expectation value on the -axis is given by the average movement per stone = run / A . The standard deviation is the same as for diffusion. In Figure 6 the result of a simulation is shown. One can clearly recognize the Gaussian distribution and the expectation value = 200 according to the above equation. It should be noted that the distribution is rather asymmetric at the very beginning because no moves in negative directions are possible. However, after a certain amount of steps it becomes symmetric. It should also be mentioned, that this approach does not allow varying the velocity of a mobile carrier gas because there is no explicit carrier gas. The transport is the result of a directional diffusion towards positive -values only.

Stationary Phase
In the next step, a stationary phase in introduced mathematically. In chromatography the stationary phase is required for the separation of substances due to the different interactions between the substances and the stationary phase. The distribution of the solute between the mobile phase and the stationary phase is described by the distribution law of Nernst: Here is the equilibrium constant, ads the concentration in the adsorbed phase, and m the concentration in the mobile phase. This law gives a linear adsorption isotherm, which is the ideal case resulting in symmetric peaks in the chromatogram. Deviations from the linearity are related to asymmetric peaks with tailing or fronting. A way to mimic the retention of a substance by a stationary phase is the introduction of a probability of a chosen molecule actually to be moved. A molecule that is strongly interacting with the stationary phase may be moved with a probability of, for example, 10% while another molecule with vanishing interaction is moved with 100% probability.
Hence, the interaction is represented by the probability to move the molecule. Once a stone is chosen in the game, one draws another random number and moves the stone depending on the outcome of the second random number. Within the simulation a second random number acc ,A is drawn for substance A from 1 to acc ,A . The stone is moved only if acc ,A = 1. This corresponds to a probability for moving of 1/ acc ,A . We can implement this in the code in Table 2(Supporting Information) by introducing just two lines between lines 7 and 8: iaccA = random (1,naccA) if iaccA = 1 then continue lines 8 to 11 Here, lines 8 to 11 become part of the conditional statement. In order to simulate the separation of two substances the algorithm needs to be doubled. A corresponding program code is provided in the Supporting Information.  Figure 7 the result of such simulation is shown. The simulated peaks are in excellent agreement with the corresponding Gaussian distribution. This chromatogram is different to that of gas chromatography because it represents the distribution of the two substances along the column ( -axis) and not the signal of the detector as function of time. On the other hand, the paper chromatography representation at the bottom actually corresponds to an experimental paper chromatogram. The retention time can be calculated from the velocity of the peak maxima and the length of the column from injection to detector at the end of the abscissa. It is obtained in reduced units and may be related to the real retention time if the length of the stationary phase is related to the real column length and the number of steps is related to the real time.

Simulation of Lighter Gas Chromatogram
The experimental realization of gas chromatography in upper secondary education is challenging but it is possible with limited resources. Choguill [31] has proposed a simple self-made chromatograph designed for high school. He used nitrogen as carrier gas, a glass column filled with firebrick soaked with heavy silicon oil, and used a thermal conductivity cell as detector. Cowan and Sugihara constructed an inexpensive apparatus for the separation of mixtures of organic substances [32]. For detection, they measured the change in temperature above the hydrogen flame by means of a thermocouple. McLean and Pauson [33] used a similar apparatus but detected the emerging chlorinated organic substancs by a copper wire that gives a green flame of volatile copper halides. Another low cost gas chromatograph suitable for high school has been developed by Kappenberg [34]. Lorke and Sommer [35] investigated how chromatography is taught in high school. They reported that chromatography is often introduced by paper chromatography in an early stage when students have not yet a concept of interactions and conclude, that it should repeated on a higher level in upper secondary education.
In order to simulate a system that we can obtain in an experiment we can use lighter gas [34]. For setting up the simulation it is necessary to relate the experimental retention time to the acceptance parameters acc ,X of component X. It turns out that a linear relation of the retention time to acc ,X is suitable. In order to avoid large values for acc ,X and hence lengthy simulation runs, all retention times can be divided by a constant.
For example, the chromatogram depicted in Figure 8 is based on the retention times 58s (ethan), 75s (propane), 107s (methylpropane), and 123s (n-butane) [34]. Here we divide these values by five and get for acc ,X the values 12, 15, 21, and 25 respectively. The experimentally obtained mole fractions of the substances are multiplied with the total number of molecules in the simulation tot = 3 • 10 4 . In Figure 8   In order to construct the chromatogram on the time axis we have to position a detector at some point along the column. In experiments, the detector is at the end of the column but if we want to observe the chromatogram on both axis it should be placed significantly before the end of the column. The number of successful moves divided by the total number of molecules in the system appears to be a useful measure for the time. Still the choice of the time axis is arbitrary and in the end it has to be related to the real time axis. If we divide the expectation value of the Gaussian distribution on the time axis (at the maximum) by the above given experimental retention time r we get a constant of about 3,63 for this specific simulation. Hence, the real time in seconds is the simulation time sim divided by 3,63. However, the peak area on the time axis does not correspond to the peak area on -axis . Since the peak area is the amount of molecules they should be identical on both axis. Here we find for the simulation shown in Figure 8 a linear relation between the ratio of the peak areas and the simulation time given by / = 0,0028 • sim (ranging from 0,6 to 1,3 for the four substances). This may be explained by the fact that a "faster" substance passes the detector faster and appears therefore slimmer on the time axis as a slower peak. This can be observed very clearly in simulations with two substances of equal amount. Here we have to be reminded that the flow of the substances in the column is caused here by directional diffusion only. There is no contribution by the velocity of a flowing carrier gas for keeping things simple. Adding a fast mobile phase would diminish this effect. Still the simulation chromatogram represents qualitatively the experimental one.

Deviation from the Linear Adsorption Isotherm
If diffusion is symmetric with respect or the center of gravity, the obtained peaks are symmetric as well. In the limit of symmetric Gaussian peaks, the single moves are stochastically independent from each other. Hence, there is no interaction between the solute molecules. However, deviations from the linear adsorption isotherm can lead to asymmetric peaks. If the peak maximum is shifted to the left on the time axis, the asymmetry is called tailing. We can model tailing by a non-linear adsorption isotherm such as: For = 0 we recover the Nernst equation 10. A negative deviation from the linear distribution law leads to an overloading of the stationary phase, i.e. not all molecules can be absorbed, and hence to faster moving of the substance in the column. To mimic this process, the fraction m of molecules in the bin at given -position is calculated. This is the number of molecules A at givenposition divided by the total number of molecules of type A. For consistency the value of m should be lower than the m -value at the maximum of this specific adsorption isotherm. So the higher the local concentration m in the mobile phase the faster the molecules are moved and hence the lower acc must be. Dividing equation 11 by m gives an effective equilibrium constant * = ( m + 1) that can be translated to the acceptance parameter: acc * = acc ( m + 1). For a negative -value this is a linearly decreasing function corresponding to increasingly faster movement of the molecules. All we need to do is to replace acc in the conditional statement in section 2.3 by the expression for acc * leading to the code line: iacc=randint (1,max (1,int(naccAN*(w*nposA[xpos[ird]]/ molA+1)))) Here the commands int and max make sure that we use an integer number not smaller than one. It is necessary to chose acc -values (naccAN) in the order of multiples of ten to get a reasonable effect. We can observe significant peak tailing for example with = −2,5 . It should be noted that tailing is obtained only after transformation to the time axis. Along the column, the peak is mirrored. The asymmetry may be numerically described by ( √ − )/ where and are calculated from the simulation data of the peak. It should be stressed that the approach for asymmetry here is just qualitatively to show the effect only.

Conclusion
Modeling and simulation has become very important in chemical science. It should be included in secondary education at a reasonable level. Students should get the opportunity to gain knowledge about the fact that computational approaches and mathematical modelling are relevant in chemistry. To be efficient, modelling should be applied to topics that are already present in the curricula in secondary education. Chromatography is potentially such a topic and can in addition be employed as a context for teaching basic concepts such as molecular interactions and diffusion. It furthermore allows linking chemical topics to mathematics in upper secondary education. It is also possible to introduce the binomial distribution in mathematics with the simulation game for diffusion. From the mathematical point of view, such modeling in science provides a justification using a real scientific context. It is therefore suggested to enrich secondary education to a reasonable extend by modeling and simulation in the best case complementary to student experiments.
For = 1/2 we get: Here can be understood as the number of steps related toone molecule. In order to calculate we have to divide the total number of steps by the number of moving molecules, which yields: