Diffusion as a Random Walk

A random walk is a sequence of discrete steps in which each step is in a random direction and has a fixed step length [1]. A perfect example is the “drunkard's walk,” where a drunk person staggers left and right with equal probability while also moving forward. Many physical processes are stochastic in nature, and random walk is the perfect model to simulate these in-silico.

One example is Brownian motion. Brownian motion of particles is a continuous random walk process due to random thermal fluctuation. Brownian motion can be viewed as each particle moving independently and randomly, where “independently” means each particle follows its own path without influence from other particles and “randomly” means each particle's position at next time step is uncorrelated to its past positions. To capture Brownian motion using random walk, each walker in the simulation represents a particle in a physical system.

Particles subjected only to Brownian motion have erratic movements at microscopic level but show net movement of descending the concentration gradient when observed on a macroscopic level. Therefore, diffusion can be seen as a macroscopic phenomenon driven by microscopic Brownian motion. A random walk with a sufficiently large number of “walkers” captures the diffusion pattern after the simulation is left to evolve for sufficient time. The connection from random walk can be drawn from the physical nature of the process, and a robust mathematically derivation can also be cleverly shown with few lines of algebra.

The movement of NiSkin from the surface of the skin to the dermis is primarily described by diffusion, but that does not mean we can boldly ignore other biological parameters. After careful consideration, we include two biological parameters: permeability coefficients and half-life of NiSkin. Considering these factors in our random walk model allows for a more realistic simulation.

Permeability of Skin

As the first line of defense against infections, the skin acts as a barrier against many microorganisms and molecules. Though TD-1, the transdermal peptide included in NiSkin, will facilitate the delivery of NiSkin deep into the dermis. However, the upper layer of the skin still poses a significant barrier to NiSkin, exerting non-negligible influence over its diffusion pattern. To quantify the skin permeability of NiSkin, we use permeability coefficients which relate solute flux to the concentration gradient across a membrane. [2] This permeability coefficient effectively defines the rate of particle exchange through the membranes. [3]

Chemical Stability of NiSkin

Another factor we consider is denaturation of NiSkin. Like all proteins, activity of NiSkin decreases over time and is reliant on factors like pH and temperature. Nisin, the antimicrobial component of NiSkin, is responsible for targeting MRSA, and its activity level is critical for the overall effectiveness of NiSkin. As NiSkin is used topically on the skin, nisin’s activity has to be evaluated within the skin’s environment. The environmental factors we consider are pH and temperature. For typical human skin:

• The skin is slightly acidic at its surface (pH 4.1-5.8), becoming more neutral deeper into the skin. [4]
• The skin has a temperature of 37°C.

Based on literature describing nisin's chemical characteristics and activity, we can use the half-life of the nisin protein to quantify the retention of activity of NiSkin as NiSkin diffuses through the skin.

Diffusion Coefficients

The diffusion coefficient \(D\) is a physical constant dependent on molecule size and other properties of the diffusing substance as well as temperature and pressure. [5] Diffusion coefficients describe the amount of a particular substance that diffuses across a unit area in 1 second under the influence of a gradient of one unit and with units of distance²/time. As there is no existing literature on our novel recombinant protein NiSkin, we approximate the diffusion coefficient of NiSkin in water based on molecular size of NiSkin.

Our NiSkin proteins, depending on the construct version, have molecular weights between 5 and 8 kDa. One protein with similar molecular weight is equine heart cytochrome C (1HRC), which has a molecular mass of 12.37 kDa. 1HRC's diffusion has been measured experimentally and has a value of \(1.3\times10^{-10}\) m²s^-1 in water at room temperature (25°C). [6] We approximate Niskin's diffusion coefficient to the same order of magnitude as 1HRC's, which is \(1\times10^{-10}\) m²s^-1 \(\sim100\) µm²s^-1.

Now we have diffusion coefficients of NiSkin in water. The next step is to approximate diffusion of NiSkin in the skin. We can safely assume that NiSkin will diffuse much slower in skin than water. One factor responsible for the slower diffusion is porosity [8], which is a measure of the void spaces in a material. We can approximate the diffusion coefficient of NiSkin using constants in the literature and a published equation relating diffusion in water to diffusion in other mediums through porosity [6]. This approximation yields the diffusion coefficient of NiSkin in stratum corneum to be 1 µm²s^-1 to the order of magnitude, and the diffusion coefficient of NiSkin in the variable epidermis and dermis to be 10 µm²s^-1.

Permeability Coefficients

Prior research has shown that the permeability coefficient \(\kappa\) is influenced by octanol-water partition coefficient \(K_\mathrm{oct}\). \(K_\mathrm{oct}\), also commonly known as LogP on a logarithmic scale, measures the ratio of a substance's concentration in the octanol phase to its concentration in the aqueous phase of a two-phase octanol/water system, which quantifies how lipophilic or hydrophilic the substance is. A partition coefficient larger than 1 means the substance is more lipophilic vs. hydrophilic, and vice versa. LogP, which is \(\mathrm{Log}(K_\mathrm{oct})\), is usually measured experimentally, but with no NiSkin sample available at the time of modeling, we decided to predict LogP of NiSkin using an existing Deep Neural Network (DNN) model. [9] The predicted LogP value for our NiSkin construct is 3.982, which means that our NiSkin is more lipophilic than hydrophilic.

We use an overall effective permeability coefficient \(\kappa\) that treats all the layers above the dermis as a whole instead of assigning each layer a distinctive value. The permeability coefficient \(\kappa\), relating to solute flux to the concentration gradient across the membrane, is expressed mathematically by equation \[\kappa=K_{m} \cdot D_{m}/\delta,\] where \(K_{m}\) is the membrane/water partition coefficient, \(D_{m}\) is the solute diffusion coefficient within the membrane, and \(\delta\) is the diffusion pathlength [2]. Hence, in our case:

• \(K_{m}\): stratum corneum/water partition coefficient, because stratum corneum is the most significant barrier along the diffusion path of NiSkin. \(K_{m}\) is calculated from our predicted LogP value of NiSkin.
• \(D_{m}\): diffusion coefficient in the epidermis, which constitutes most of the thickness above the dermis.
• \(\delta\): effective diffusion path when molecules path through stratum and epidermis.

Sanity Check: Does this approximation make sense?

We include a controlled group, human epidermal growth factor (hEGF) and human epidermal growth factor with TD-1 attached (TD-1+hEGF). Research has shown that by attaching TD-1 to hEGF, the transdermal efficiency can increase by 20-fold. [11] Hence, our results should be consistent with this experimental finding, which means TD-1+hEGF should have a larger predicted permeability coefficient than hEGF alone.

We also perform the calculation on nisin without the transdermal peptides TD-1 attached. The permeability coefficient should be lower than NiSkin because TD-1 facilitates transdermal diffusion.

Our results are listed below and the results agree with existing experimental data and our hypothesis that TD-1 will improve nisin transdermal ability.

Half-life of NiSkin

Assuming denaturation of NiSkin is random, which means that each protein molecule has equal probability of denaturing at each unit time, the denature process can be modeled like a radioactive decay process. The number of radioactive particles \(N\) undergoes exponential decay, which is described by the equation \begin{equation} N(t)=N_{0}e^{-\lambda t}\tag{1} \end{equation} where \(N_{0}\) is initial number of active particles, \(\lambda\) is often known as the decay constant, and \(t\) is the time. Research has shown that under neutral (pH = 7) environment, nisin retains about 20% of activity after incubation at 180 RPM for 3 hours at 37°C [12]. Using Eqn. (1), we can easily calculate the half-life of nisin, which is the time it takes for nisin to reduce to 50% activity.

Assumptions

British statistician George Box once said, “all models are wrong, but some are useful.” All models are wrong because models are overly simplified descriptions of the real world, so there are limits to their accuracy and applicability. Useful models, however, manage to capture the essence of what is happening in the real world by making careful assumptions that simplify the problems yet provide useful information.

The assumptions we make to keep our problem simple yet capture the overall diffusion pattern:

• We assume the layers above dermis have a single collective permeability coefficient.
• We assume the denaturation of nisin follows exponential decay.
• We assume that NiSkin diffuses at equal speed in all directions.
• We assume the clearance rate of NiSkin due to blood flow at dermis is negligible, where clearance rate means the the rate that NiSkin disappears from the system.
• We assume NiSkin does not interact with molecules present in the skin.
• We assume the surface of the skin as an impermeable membrane, i.e., NiSkin cannot diffuse into the air.
• We assume that the bacteria causing the infection are halfway deep into the dermis.

Random Walk Parameters

Let's now start a random walk simulation. The definition of random walk is a sequence of discrete steps in which each step is in a random direction and has a fixed step length. Hence, two factors are important here, sequence of steps, which is determined by time, and fixed step length.

• Let's assume we have a time series \(t=(t_1,t_2,\ldots,t_n)\), where all the timesteps are evenly spaced. So, the same time of ∆𝑡 has passed between \(t_1\) and \(t_2\), and between \(t_2\) and \(t_3\), and so on.
• Let's assume that each walker travels the fixed step length of \(\Delta x\).

All the workers in the simulation follows these exact same steps:

1. At time = 0, the walker is put into the simulation and ready to start from its starting coordinate.
2. At time = \(t_1\), a new random position that is \(\Delta x\) away from walker's original position is proposed.
3. At time = \(t_1\), the walker accepts the proposal and moves to the new position.
4. At time = \(t_2\), a new random position that is \(\Delta x\) away from walker's original position is proposed.
5. At time = \(t_3\), the walker accepts the proposal and moves to the new position.
6. At time = \(t_3\), a new random position that is \(\Delta x\) away from walker's original position is proposed.
7. ...

These steps are repeated as many times as needed to let the system evolve with time.

The mathematical derivation proves that the relationship between random walk and diffusion equation is non-trivial, but you may be wondering how step size \(\Delta x\) and time \(\Delta t\) are related to diffusion coefficients. It turns out that these three quantities are nicely related by one equation, \[\Delta x=\sqrt{6D\Delta t}\] which can also be proven mathematically.

This random walk simulation belongs to a broader class of computational algorithm, Monte Carlo methods. Monte Carlo methods are computational algorithms that generate suitable random numbers to obtain numerical results and the random number in the case of random walk is the new position the walker will move to.

Different versions of Monte Carlo have been developed to address different problems. To incorporate both half-life of NiSkin and the permeability coefficient into random walk, we choose to use Metropolis Monte Carlo, where each walker has a probability of accepting the proposed next position rather than always moving to the proposed position. How is this probability implemented in the random walk? For example, let's say the probability we want to enforce is 0.7. A random number between 0 and 1 is generated, and if the value is smaller than 0.7, which happens 70% of the time, the walker accepts the proposal and moves to the proposed position. Otherwise, it will stay in the same position until the next step, which happens 30% of the time. Using some mathematical tricks and equations proven in literature, we include half-life of NiSkin and permeability coefficient as two probabilities in the Metropolis Monte Carlo simulation.

The 3D random walkers in the simulation represent billions of NiSkin molecules in real life. To simulate NiSkin applied onto skin in an area, we map random walkers that originate from one point to a square of 1 cm² We divide each side by \(n\) times, and now we have \(n^2\) grids inside the square. Assuming \(n\) is very large, all the walkers starting at every vertex of the grids can be seen as covering the entire 1 cm² area.

Figure 1. The diagram of modeling setup.

The simulation aims to answer how the concentration of NiSkin changes inside the skin, and hence, we need to monitor all the trajectories of the walkers to keep tracks of the position of NiSkin at different time points.