Our engineered protein NiSkin is intended to be used as a topical treatment to effectively treat cellulitis. Due to time constraints and ethical concerns, not all of the questions pertaining to NiSkin effectiveness can be answered by consulting experts or experimenting during wet lab. This list of questions includes:

- How much NiSkin penetrates the stratum corneum?
- What is the concentration of NiSkin in the dermis after applying the cream?
- How often do we need to apply NiSkin so that the concentration remains above minimum inhibition concentration (MIC)?

We turn to modeling to tackle these questions. With reasonable assumptions and careful selection of parameters, *in-silico* modeling can simulate reality with reasonable accuracy.

Imagine applying NiSkin onto your forearm. How will the NiSkin move to the cellulitis infection in the dermis? The driving physical process responsible for this movement is **diffusion**. Just like how a food's smell travels across the room as particles move from an area of high concentration to an area of low concentration, NiSkin diffuses from the surface of the skin to the dermis down a concentration gradient.

The question lies in how to capture this diffusion process with modeling.

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.

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]

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.

After selecting the parameters to include in our models, the next step is to determine the values of these parameters, either directly from existing literature or by making reasonable approximations.

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^{2}/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^{2}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^{2}s^{-1} \(\sim100\) µm^{2}s^{-1}.

Other proteins around the same molecular weight also have experimentally measured diffusion coefficient in water at the same order of magnitude. For example, egg white lysozyme (1BWI) has a molecular weight of 14.33 kDa and a diffusion coefficient of \(1.12 \times10^{-10}\) m^{2}s^{-1} [6].

Another way to estimate diffusion coefficient is using Stokes-Einstein's equation. By approximating NiSkin to a perfect sphere with a density of 1.35 g/cm^{3} [7], the diffusion coefficient calculated has the same order of magnitude.

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 ^{2}s^{-1}** to the order of magnitude, and the

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.

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.

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.

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.

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:

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

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^{2} 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^{2} area.

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.

The code for random walk simulation is written in C and executed on Rivanna, the high-performance computing cluster (HPC) at University of Virginia. \(n\) is taken to be 100 for a small enough grid, but still stay computationally feasible. 10,000,000 walkers are simulated in the system, with each walker representing NiSkin molecules in the real world. To capture how the concentration of NiSkin changes with time, the random walk simulation is left to evolve for 8 hours.

The number of walkers at each layer of the skin is counted every 10 minutes. The results are shown in Figure 2. We can see that the number of walkers entering the dermis region in the simulation peaks at around 1 hour.

We also modeled different reapplication intervals of NiSkin and compared the distribution profile. Reapplication intervals of 1, 2, 3, and 4 hours are shown below (Figure 3). The distribution profile given by 1-hour and 2-hour interval shows that the number of walkers increase for all layers, which means the concentration within each layer will build up over time, reaching dosage that could cause allergic reaction. The reapplication with 4-hour intervals gives a more regular pattern, that suggests the concentration at each layer will be better under control.

When calculating LogP value and the permeability coefficient, we show that our NiSkin construct has much higher permeability coefficients than nisin without transdermal peptide TD-1. This result suggests that our NiSkin construct is highly like to demonstrate transdermal property.

From random walk simulation, we can confirm that around 20% of walkers enter the dermis after 1 h, which again shows that our NiSkin construct is very feasible. Our modeling also shows that reapplication after 4 hours can bring the number of walkers in the dermis region to approximately the same concentration before, without the chance of overdosing.

To improve our model, we can incorporate additional parameters to increase the realism of our model. We use models to find the concentration at a certain depth of the skin to gain more insight into the diffusion process. For example, how concentrated our NiSkin cream needs to be to allow the concentration at the dermis to reach the minimum inhibition of the pathogen in dermis. We on one hand need to make sure the concentration of NiSkin is high enough to reach the optimal dosage. On the other hand, we cannot increase concentration indefinitely, because research has shown that nisin could be an irritator or sensitizer at high concentrations, which refers to a substance that produces allergic reaction after repeated contact with skin. Future work will focus on finding the optimal concentration that satisfied the concentration limit research has specified.

- Random Walk—From Wolfram MathWorld. (n.d.). Retrieved October 12, 2023, from https://mathworld.wolfram.com/RandomWalk.html
- Potts, R. O., & Guy, R. H. (1992). Predicting skin permeability. Pharmaceutical Research, 9(5), 663–669. https://doi.org/10.1023/a:1015810312465
- Alemany, I., Rose, J. N., Garnier-Brun, J., Scott, A. D., & Doorly, D. J. (2022). Random walk diffusion simulations in semi-permeable layered media with varying diffusivity. Scientific Reports, 12(1), Article 1. https://doi.org/10.1038/s41598-022-14541-y
- Proksch, E. (2018). pH in nature, humans and skin. The Journal of Dermatology, 45(9), 1044–1052. https://doi.org/10.1111/1346-8138.14489
- DIFFUSION COEFFICIENT. (n.d.). Retrieved October 12, 2023, from https://www.thermopedia.com/content/696/
- Wang, J., & Hou, T. (2011). Application of Molecular Dynamics Simulations in Molecular Property Prediction II: Diffusion Coefficient. Journal of Computational Chemistry, 32(16), 3505–3519. https://doi.org/10.1002/jcc.21939
- Densities of molecules, proteins, organelles,—Various—BNID 101502. (n.d.). Retrieved October 12, 2023, from https://bionumbers.hms.harvard.edu/bionumber.aspx?id=101502
- Liu, J., Ding, W., Ruan, R., Zou, L., Chen, M., Wei, P., & Wen, L. (2017). A Theoretical Study on Inhibition of Melanoma with Controlled and Targeted Delivery of siRNA via Skin Using SPACE-EGF. Annals of Biomedical Engineering, 45(6), 1407–1419. https://doi.org/10.1007/s10439-017-1825-5
- Chen, Yan-Kai, Steven Shave, and Manfred Auer. 2021. "MRlogP: Transfer Learning Enables Accurate logP Prediction Using Small Experimental Training Datasets" Processes 9, no. 11: 2029.
- Adhikari, M. D., Das, G., & Ramesh, A. (2012). Retention of nisin activity at elevated pH in an organic acid complex and gold nanoparticle composite. Chemical Communications, 48(71), 8928. https://doi.org/10.1039/c2cc34653b
- Alemany, I., Rose, J. N., Garnier-Brun, J., Scott, A. D., & Doorly, D. J. (2022). Random walk diffusion simulations in semi-permeable layered media with varying diffusivity. Scientific Reports, 12(1), Article 1. https://doi.org/10.1038/s41598-022-14541-y