Introduction

A model's purpose is to help us understand a phenomenon, predict its behavior and guide future action. This model serves to provide an understanding of the basis of our project, the Phage Litmus technology. It also showcases a journey from the bacteriophage nanostructures, to the Volatile Organic Compounds they interact with, to what needed and what still needs to be investigated in order for this technology to be applied for our diagnostic purposes.

As stated in the Project Description, the shift in the structural coloration of the Phage Litmus upon interaction with a ligand (in our case, a Volatile Organic Compound) is caused by M13 phage bundle swelling. This swelling can be quantified as the difference in bundle diameter (ΔD) and/or the difference in phage matrix thickness (ΔT) (1).

We will assume that each ligand molecule interacting with the receptor (in our case, pVIII) contributes equally to the swelling effect. Furthermore, each wafer will contain an extremely large number of bacteriophages (all exhibiting the same pVIII modification), which, combined with the large amount of pVIII copies per bacteriophage (~2700), suggests a measurement procedure with high resolution, capturing a broad spectrum of states from zero to saturation. These statements are supported by experimental data done on the detection of explosive substances in the air (2) and relative humidity (3), suggest a linear or close to linear correlation between ligand amount and phage swelling, up to the point of receptor saturation:

\[n_{\text{ligand}} = \kappa \cdot \Delta T \quad [1] \] \[n_{\text{ligand}} = \lambda \cdot \Delta D \quad [2]\]

Where κ, λ are arbitrary scalars ∈ R⁺.

Our testing will take place in an enclosed space of fixed volume.

\[ n = C * V \]

This means we can substitute moles with concentration for equations [1],[2]

\[C_{\text{ligand}} \cdot V = \kappa \cdot \Delta T \] \[ \Rightarrow C_{\text{ligand}} = \frac{\kappa}{V} \Delta T \quad [3] \]

\[C_{\text{ligand}} \cdot V = \lambda \cdot \Delta D \] \[ \Rightarrow C_{\text{ligand}} = \frac{\lambda}{V} \Delta D \quad [4] \]

Since V is fixed, ^κ⁄_V and ^λ⁄_V are scalars, thus C is proportional to T and D. Structural color shift is a direct result of phage bundle swelling. If we consider only one ligand interacting with our phages at a time, and accounting for equations [3] and [4], we can assume that, in those ideal conditions, phage bundle swelling (expressed in ΔΤ and ΔD) and color difference (ΔColor for now) exhibit a correlation.

\[f(\Delta T, \Delta D) = \Delta \text{Color} \quad [5]\]

One should note that this ΔColor is an abstract concept for now, and highly dependent on the approach and model chosen.

By combining [3],[4],[5], we can define:

\[f\left(\frac{V}{\kappa} \cdot C_{\text{ligand}}, \frac{V}{\lambda} \cdot C_{\text{ligand}}\right) \]

\[ = \Delta \text{Color} \quad [6] \]

Which means that there exists a potentially discoverable correlation between the concentration of our ligand and the colorimetric change. Determining the exact function is a subject for future research, which would be possible only after constructing the device and acquiring sufficient real-world data by testing it on many different ligand concentrations.

Given that most color difference definitions involve the measurement of distances within a color space, one method for quantifying them is the utilization of the Euclidean distance metric. When one possesses an RGB (red, green, blue) tuple and seeks to computationally determine the color difference, a straightforward approach involves considering the R, G, and B values as linear vectors that define the orthogonal color space system, which in our case is three-dimensional.

This is not a true space in the mathematical sense, since it is bounded by the constraints of the color model (in our case, 0 to 255 for each axis).

The difference in coloration (Color Distance, Cd) can be geometrically defined as follows,

\[Cd = \sqrt{(R_2 - R_1)^2 + (G_2 - G_1)^2 + (B_2 - B_1)^2} \quad [7]\]

where \( R_1,G_1,B_1 \) and \( R_2,G_2,B_2 \) are intensities of red, green and blue before and after the ligand-receptor interaction, respectively.

This distance is a mathematical representation of the ΔColor we mentioned earlier. It should be noted that the Euclidean distance metric is far from the ideal approach, since it forces a multidimensional input into a one-dimensional positive number line, while also containing many possible instances of different inputs producing identical outputs. A lot of information is therefore lost, and two distinct conditions may fail to be differentiated.

To minimize this information loss, a better approach would include calculating the difference of the 2 RGB vectors, creating a 3-dimensional vector output.

\[f\left(\frac{V}{\kappa} \cdot C_{\text{ligand}}, \frac{V}{\lambda} \cdot C_{\text{ligand}}\right) \]

\[= \vec{(R_1, G_1, B_1)} - \vec{(R_2, G_2, B_2)}\quad [8] \]

We should note that the RGB color model is not integral to this process: Any other color model that is derived from distinct axes, as numerous as these axes could be (e.g. 4 in CMYK), can be incorporated in the above function definition.

Since our output will be captured with a smartphone camera and converted into RGB, or any other color model values, some kind of quantization will take place. However, since \(256\) integer values for each color create \(256^3=16777216 \) possible colors, we can substitute the continuous shift in structural coloration with the RGB-converted colorimetric output with only limited information loss.

Once sufficient testing data is acquired and the above function (equation [6]) is better understood, substituting the RGB values will enable us to reliably semi-quantify the concentration of our target compound in the sample, thus giving us valuable insight into the target’s propensity for developing Parkinson’s Disease. However, due to the complexity of the optical phenomena and the sample composition, a reliable form of semi-quantification may never be possible. Therefοre, more indirect methods of analysis could be investigated.

In order to differentiate the colorimetric results between PD patients and healthy participants it is important to study the color boundaries of healthy and PD results. Thus, other color channels/coordinate systems should be examined as well, like YIQ or the aforementioned CMYK. The conversion formulas across different channels will facilitate this boundary colorimetric result study.

As we mentioned above, there hasn’t been a colorimetric Phage Litmus result analysis for PD-related VOCs. So, the colorimetric boundaries might not be well established by any of the traditional color channels and linear or non-linear transformations might be imperative to be applied on the color channels/coordinate systems. This fact, coupled with the fact that our colorimetric results consist of multiple litmuses of possibly independent color gradients, make it really hard for the boundary to be simple enough to be manually calculated or for the necessary linear or nonlinear dimensionality transformation functions to be tracked.

For the reasons listed above, a Deep Learning pipeline could be implemented, to uncover those hidden transformations and boundaries of PD patients and healthy controls. The deep learning model goals are tailored to the boundary and dimensionality problems we encountered, because they are comprised of many layers of nonlinear functions that aim to make the samples separable linearly.

At this point, after having ventured this far, we should remember the goal of our project: Creating a sensor for detecting specific compounds. This sensor is colorimetric, so we aim for the maximum color change when detecting our compounds of interest, and the minimum change when detecting anything else. For the purposes of standardization, we would appreciate the behavior of the colorimetric response function to be as clear as possible. Unfortunately, the sample of evaporating substances from skin sebum is mixed, with most compounds (probably) not related to Parkinson’s Disease. Furthermore, our sensor cannot be perfectly specific for its target compound. These limitations would create complex, unpredictable responses and noise in the phage’s colorimetric output according to the ligand’s concentration levels, as other compounds of similar structure would interfere with our biomarkers binding to their respective receptors.

To eliminate as much of this noisy output, there are certain goals we must strive towards:

• We must maximize the affinity of our receptors for their corresponding ligands.
• We must minimize the affinity of our receptors for any other ligand in the sample.

Therefore, we needed to find receptors that would exhibit high affinity for the Volatile Organic Compounds of interest. In addition, we needed to find a method for filtering out receptors with high affinity for unrelated compounds contained in human sebum. As described in the literature (1,2,3), the pVIII protein part responsible for interacting with the Volatile Organic Compounds is the area proximal to the N-terminus, whose free component (not participating in the α-helix) is composed of 5 amino acids.

Literature research was conducted on proteins known to interact with our ligands of choice to identify an appropriate modification. Databases of Odorant Binding Proteins (5) were investigated, but they did not provide us with an amino acid ligand-binding sequence specific for our compounds. Thus, we had to find it ourselves.

Phage display

Phage litmus techniques have always been reliant on phage display to determine a consensus ligand-binding sequence (6). This is a wet lab procedure, but describing it in brief will provide some reasoning behind our choices further below.

Phage display technology is an in vitro screening technique for identifying ligands for proteins and other macromolecules. At the crux of phage display technology is the ability to express peptide or protein sequences as fusions to the coat proteins of a bacteriophage.

A typical phage display procedure is conducted as follows:

Target ligands are first immobilized to the wells of a microtiter plate. Many genetic sequences are expressed in a bacteriophage library in the form of fusions with the bacteriophage coat protein, so that they are displayed on the surface of the viral particle. The protein displayed corresponds to the genetic sequence within the phage.

This phage-display library is then added to the dish and after allowing the phage time to bind, the dish is washed. Phage-displaying proteins that interact with the target molecules remain attached to the dish, while all others are washed away. (7)

Due to limited resources, time and experience, we decided to conduct our biopanning by using computational methods. To do that, we had to simulate the procedure.

We mentioned our need for a process that would determine the binding affinity of a ligand against a receptor. Upon consulting with several bioinformaticians (listed in the Attributions and Integrated Human Practices pages), we concluded that Molecular Docking software should be used for this purpose.

Molecular Docking

Molecular docking is the study of how two or more molecular structures fit together. In a simple definition, docking is a molecular modeling technique that is used to predict how a protein (receptor) interacts with small molecules (ligands). Several software packages exist for such purposes.

Upon further discussion and research, we came upon Autodock Vina (8). Compared with other molecular docking software, Vina provided us with a fast and accurate method for screening through a large combinatorial space of possible receptors.

In trying to emulate phage display, the following pipeline came to mind:

• We will create many different amino acid combinations (oligopeptides).
• We will insert the said combinations proximal to the N-terminus of the pVIII protein
• We will perform molecular docking using the modified protein as a receptor and each of the four Volatile Organic Compounds of interest as ligands.
• Proteins with high affinity results will be chosen for further analysis, and the others will be discarded.

Initially, to make our simulation as realistic as possible, we tried to test the whole modified pVIII protein against our ligands. This would include inserting a random oligopeptide into the protein’s sequence, folding it (using AlphaFold) (9,10) and then docking.

This posed several challenges. First of all, folding is not a fast procedure. Even with AlphaFold2, which accelerated the process considerably, it would take minutes per protein. To discover a sensitive and specific protein modification, a very large number of receptors would need to be tested, and time was of the essence. Furthermore, we would face search grid-related issues. In molecular docking software, including Autodock Vina, users need to specify the dimensions of the space the algorithm would be asked to search. Loading a large protein would position the N-terminus in an unpredictable manner, risking the generation of results from parts of the protein not exposed to the outer surface of the bacteriophage (useless). Thus, another solution was found.

To simplify and accelerate the procedure, we decided to test small peptides against our ligands. This addressed both problems: We would avoid the need for folding, as oligopeptides do not usually adopt unpredictable 3D conformations, only needing preparation and 3D structure optimization which is relatively fast and easily automated. We would also considerably simplify the generation of grids, as a relatively small one could cover the whole peptide-ligand complex if positioned in the center (0,0,0). It was also considerably easier to write code for automated grid placement, which would find the center of mass and extend the grid to the end of the peptide. This code, along with every other function we used and wrote to perform and automate the docking procedure is provided in our Software page and GitLab.

So, we decided to start small: tripeptides. All possible tripeptides (8000) were generated, prepared, and converted to PDBQT format (Vina standard). During preparation, water molecules are removed, polar hydrogens are added, and partial charges are assigned. Removal of water molecules is somewhat of a controversial topic in molecular docking. We decided to do it, since it is standard practice, and because our sensor is designed to work in atmospheric conditions and not in solution. Ligand files were also downloaded from PubChem and prepared accordingly.

Then, the tripeptides were tested against each ligand in a loop. Hippuric acid and perillaldehyde took under 10 seconds per peptide, while eicosane and octadecanal took around 30 seconds.

Once this process was completed, binding affinity results were organized and analyzed. Python code was written to compare the docking results between different ligands. Since perillaldehyde sebum amounts are lower in Parkinson’s Disease, in contrast with the other three compounds (11), it was of crucial importance to choose a receptor that would not bind with the other three VOCs to a considerable extent. For semi-quantification purposes, specificity between the three PD-positive biomarkers was important, but not critical.

We mentioned a while back that it was important to establish specificity for our receptor when it came to other VOC present in human sebum. To achieve this, we delved into the literature and found several papers documenting them. We chose 21 compounds found on the upper back, which is also the body region from which we plan to obtain our samples (12). The reason we chose this body region is that it matches the sampling region of the original paper on PD-related sebum VOC (11). As we did with the other ligands, we downloaded and prepared those compound files as well. Since checking those VOC for every possible receptor would multiply our docking time by 21, we chose to do this only for our top contenders as a final evaluation.

These tripeptide results were not satisfactory as of both our selection criteria (in essence, high specificity and sensitivity). To extend our combinatorial space further in the hopes of discovering better candidate sequences, we increased our peptide length to 4. Thus, the number of peptides to be created and tested increased from \(20^3 = 8000\) to \(20^4= 160000\). This posed certain practical issues. We were not able to secure a cluster for our computations, relying instead on our home computers. The docking times mentioned above resulted in an estimated duration counted in months. This was prohibitive, as orders and wet lab experiments would need to be initiated in a timely manner. This forced us to focus on certain ligands (perillaldehyde and hippuric acid) more so than on others (octadecanal and eicosane). We did not abandon them, but we tested fewer peptides against those ligands.

In addition, we recognized the need for targeted peptide testing. When looking through our more promising results, we noticed certain amino acids or amino acid patterns being overrepresented (eg. tryptophan in eicosane results), and directed our efforts accordingly.

We also noticed patterns of hydropathicity among our top contenders. They tended to tilt towards the side of hydrophobicity, especially when it came to eicosane and octadecanal, perhaps not surprisingly. A plot of all peptides we analyzed ranked by affinity and with the Grand Average of Hydropathicity (13) calculated is shown on the below:

Our needs for timely results with limited computational strength, along with the patterns we noticed in the data provided inspiration for a targeted docking algorithm that would assist low-resource researchers in finding accurate results sooner when exploring large receptor spaces. Further information about our efforts in developing this algorithm is presented on the Software page of our Wiki.

Wet lab protocols (14) provided an additional limitation: our chosen sequences would need to begin with A , V , D, E , G. The details surrounding this choice are displayed on the Design page of our Wiki. Results already obtained and satisfying this condition were not found to be of acceptable quality. Therefore, we created pentapeptides based on this new data.

At this point, we decided to once again use AlphaFold2 to check whether peptide results and engineered pVIII results were close enough. For this, we chose a then- candidate hippuric acid sequence (DYYAW), inserted it into the pVIII protein at the appropriate position and performed docking. It should be noted that protein modifications were guided by precisely following the wet lab protocol using software such as SnapGene, so as to emulate the actual situation as accurately as possible. This is when we noticed a discrepancy:

Folded protein results were much worse than on the peptide alone. We attributed this discrepancy to the presence of Alanine on position 1, right before the point of insertion. To validate our results, we performed docking with the hexapeptide ADYYAW and found similar results with the folded protein. This led us to the conclusion that we should test hexapeptides, starting with Alanine, followed by one of A, V, D, E, G, followed by 4 random amino acids. Of course, all possible hexapeptides (64000000) is a completely prohibitive search space for a team with our time and resource constraints. The fixed Alanine in position one limits our choices to 3200000, and the 5 choices for the second amino acid bring us down to 800000, still an extremely high number. To determine whether we should create hexapeptides based on our best already tested tetrapeptides, we compared the affinity results between tetrapeptides and hexapeptides whose last 4 amino acids matched. We found no significant correlation, so we created them randomly.

From beginning to end, we managed to test almost 300000 peptides against hippuric acid, almost 190000 peptides against perillaldehyde, over 46000 peptides against octadecanal and over 45000 peptides against eicosane, amounting to a total of almost 5 million seconds of run time, (57 days) on our main computer. In the end, the wet lab procedure did not require the limitations imposed on the second amino acid in the sequence, since another method was used (see Engineering). This meant that peptides starting with A, followed by A,V,D,E,G were severely overrepresented. A database containing our docking results is presented and described on our Results page.

We also showcase our chosen peptide sequences for each ligand of interest, which was each compared against the 21 sebum VOC and with the other three PD-related VOC, as well as the unmodified pVIII sequence. As a final test before determining our final sequences, we used AlphaFold2 one last time, to compare results between hexapeptides and folded pVIII protein. Results did not match completely, but all the previously established conditions were still satisfied.

Further Analysis

The aforementioned criteria of sensitivity and specificity were being enforced manually. To further automate the procedure, a scoring function will be developed and implemented into code, which would select candidate sequences based on the maximum binding affinity difference from other compounds in the sample, while also selecting for high binding affinity to our target compound.

To further validate our docking results, DFT and molecular dynamics were considered, so as to also dynamically evaluate our system’s (ligand-receptor complex) behavior taking into account the passage of time. DFT calculations using ORCA (15) and the B3LYP level of analysis were initiated, but not expected to conclude and be thoroughly analyzed before competition deadlines.

Now that we have investigated and modeled the behaviour of our kit once it comes into contact with the Volatile Compounds of our choice, lets take a step back. How did those compounds reach our bacteriophage-based device to begin with? In order to elucidate this process and garner invaluable information for designing and implementing our potential device, another model was developed.

This will complete the cycle and result in a modeling work encompassing our whole design, from beginning to end.

Evaporation Diffusion

Introduction

In this model, we examine the amount of time required for every Volatile Organic Compound of interest to travel from the patient’s skin sebum to the wafers in which our genetically engineered bacteriophages remain encrusted. The wafers are placed in our 3D printed kit, so that each wafer corresponds to a differently engineered phage. The kit includes 5 wafer positions in total: one for the control bacteriophages - which have not been genetically engineered -, and four positions for four different genetically engineered phages. Each of these phages will detect a different VOC.

Our kit’s structure is presented below. It has a length of 5.5 cm, a width of 3 cm, and a height of 2 cm. According to its use, our kit is placed airtight on the patient’s upper back, an area of the human body that is known to be sebaceous. The area of the human skin that shall be covered is equal to

\( \mathbf{A =length * width } \)

\[ \mathbf{= 5.5 * 3cm^2 } \]

\[ \mathbf{= 16.5 * 10^{^{-4}}m^2 }\]

The process has been divided into two parts, one subsequent of the other: evaporation of VOC’s from the skin sebum, and diffusion of the latter inside our kit.

2. Evaporation

2.1 Constructing the General Formula

Our main goal is to calculate the evaporation rate of every VOC of our interest. Our results follow Nielsen’s et al equation on the evaporation rate of a pure liquid, which was based on Mass Transfer Theory:

\( \mathbf{E = K_{g}*(C_{ii}-C_{i})} \)

\[ \begin{array}{ll} E : & \text{Evaporation Rate} \\ K_{g} : & \text{Gas phase Mass Transfer} \\ & \text{Coefficient}\\ C_{ii} : & \text{Air concentration of VOC i} \\ & \text{in equilibrium with its liquid form} \\ C_{i} : & \text{Air concentration of VOC i } \\ & \text{in the room} \end{array} \]

The concentration of each VOC in the room (air phase), before we use the kit, is very small, especially compared to the concentration of its liquid form in the stratum corneum. This means that by the time we cover the desired skin area with our kit,

\[ \frac{C_{i}}{C_{ii}} \rightarrow 0^{+} \]

so that

\[ k_{g}(C_{ii}-C_{i}) \cong k_{g}C_{ii} \quad \left [ 1 \right ] \]

One more approximation using the Ideal Gas Law yields:

\[ E = k_{g}\frac{n_{ii}}{V} = k_{g}\frac{p_{ii}}{RT} \quad \left [ 2 \right ] \]

The similar equation in correspondence with the liquid phase mass transfer coefficient is the following:

\[ E = k_{\text{l}}d*\frac{1000}{MW} \quad \left [ 3 \right ] \]

where \(d\) represents the density of the VOC in liquid form and \(MW\) its Molecular Weight.

\( \begin{align*} \text{Eicosane} & : 282.5 \, \text{g/mol} \\ \text{Perillaldehyde} & : 150.2 \, \text{g/mol} \\ \text{Hippuric Acid} & : 179.2 \, \text{g/mol} \\ \text{Octadecanal} & : 268.5 \, \text{g/mol} \\ \end{align*} \)

Table 1: Molecular Weights of VOCs of interest

\[ \begin{align*} \text{Eicosane} & : 0.789 \, \text{g/cm}^3 \, \text{at} \, 293 \, \text{K} \\ \text{Perillaldehyde} & : 0.965 \, \text{g/ml} \, \text{at} \, 298 \, \text{K} \\ \text{Hippuric Acid} & : 1.308 \, \text{g/cm}^3 \, \text{at} \, 293 \, \text{K} \\ \text{Octadecanal} & : 0.777 \, \text{g/ml} \, \text{at} \, 298 \, \text{K} \\ \end{align*} \]

Table 2: Density of VOCs of interest in liquid phase

As it will be mentioned below, the temperature at the stratum corneum is approximately 37 degrees Celsius, which corresponds to 310 K. Thus, we can calculate the density of each desired VOC in the stratum corneum, by using the following formula, which describes the relationship between the density of a VOC and the temperature:

\[ \mathbf{d_2 = \frac{d_1}{1 + \beta \left ( T_{2} - T_{1} \right )}} \]

where β is the coefficient of thermal expansion and is specific to every VOC.

For example, \[ \beta_{\text{eicosane}} = 9 \times 10^{-4} \, \text{K}^{-1} \]

It is worth noting that the liquid form of each VOC does not exist in the skin sebum as a continuous film. On the other hand, it is dissolved in the skin, and as a result we transform the eq [3] into

\[ \mathbf{E = k_{\text{l}}d\frac{C}{C_{\text{sat}}}} \]

where C is the concentration of the VOC in the uppermost layer of the stratum corneum and \( C_{sat} \) is its solubility in this layer. However, since not all of these vales are recorded in the literature, we stick with the previous equations. We combine the equations [2] and [3] to construct a formula for \( \mathbf{k_{\text{l}}} \) :

\[ \mathbf{k_{g}\frac{p_{ii}}{RT}=k_{\text{l}}d*\frac{1000}{MW}} \]

\[ \mathbf{ \Rightarrow k_{\text{l}}=\frac{k_{g}p_{ii}*MW}{1000dRT} \left [ 5 \right ]} \]

It is immediate that in order to calculate the Evaporation Rate of each VOC, one shall estimate the gas phase mass transfer coefficient \(k_{g}\). This parameter can be calculated using the United States' Environmental Protection Agency estimation:

\[ \mathbf{k_{g}=0.01756 \cdot MW^{-\frac{1}{3}} \cdot u^{\gamma} \left [ 6 \right ]} \]

Thus, substituting \(k_{g}\), we conclude that

\[ \mathbf{k_{l} = \frac{1756 \cdot 10^{-5} \cdot MW^{\frac{2}{3}} \cdot p_{ii} \cdot u^{\gamma }}{10^{3} \cdot dRT} = } \]

\[ \mathbf{ 1756 \cdot 10^{-8} \cdot p_{ii} \cdot u^{\gamma } \frac{\sqrt[3]{MW^2}}{dRT} \left [ 7 \right] } \]

Here, the exponent gamma is associated with airflow velocity in the evaporation model, and is set to be 0.5, a number that implies laminar air flow conditions.

The primary source of air movement in a closed room without mechanical ventilation is natural convection, driven by temperature differences. The speed of natural convection currents can vary depending on the temperature difference between the top and bottom of the room, but is generally quite low, often less than \(0.1 m/s \). In addition to that, there may be minor sources of air movement, which can slightly disrupt the still air. However, these disturbances typically result in very low air velocities, often less than $0.2m/s$. Thus, we solve our model with the values of u being in range of \( \mathbf{\left [ 0.1, 0.2 \right ]} \).

Having taken everything into account, we reach at the final formula by substituting \( k_{l} \). This formula will be used in the next stage of the model (diffusion):

\[ \mathbf{E = 1756 \times 10^{-5} \times p_{ii} \times } \]

\[ \mathbf{ u^{\gamma} \times \frac{MW^{-\frac{1}{3}}}{RT}} \left [ 8 \right ] \]

2.2 Evaluating the Vapor Pressure for every VOC

The temperature in the stratum corneum, which is the outermost layer of the skin, is approximately the same as the core body temperature. The core body temperature in humans is typically around 98.6 degrees Fahrenheit (37 degrees Celsius). The skin temperature can vary depending on factors like ambient temperature, blood flow, and activity level, but generally remains close to the core body temperature to help maintain homeostasis and protect the body from temperature extremes. The stratum corneum itself does not generate heat; instead, it acts as a barrier to help regulate heat loss from the body.

We can safely say that the evaporation process takes place at approximately 37 degrees Celsius. We quote the Vapor Pressure of some VOCs of our interest at 25 degrees Celsius (298 Kelvin):

Eicosane : \(4.62 \times 10^{-6}\) mm Hg

Perillaldehyde : \(434 \times 10^{-4}\) mm Hg

Hippuric acid : \(76 \times 10^{-8}\) mm Hg

Octadecanal : \(3.41 \times 10^{-4}\) mm Hg

Vapor Pressure at standard conditions (298 K)

We quote as well their Enthalpy of Evaporation, corresponding to the amount of heat energy required to vaporize a given amount of each VOC at standard conditions:

\begin{align*} \text{Eicosane} & : 102 \pm 2 \, \text{kJ/mol} \\ \text{Perillaldehyde} & : 47.5 \pm 3 \, \text{kJ/mol} \\ \text{Hippuric acid} & : 76.4 \pm 3 \, \text{kJ/mol} \\ \text{Octadecanal} & : 62.5 \, \text{kJ/mol} \\ \end{align*}

Heat (Enthalpy) of Evaporation at standard conditions (298 K)

Generally speaking, in order to find the vapor pressure of a VOC for a specified amount of temperature, one may apply Antoine's equation, which relates vapor pressure to temperature:

\[ \mathbf{\log_{10}P = A - \frac{B}{T+\Gamma }} \]

where \(A, B \) and \( \Gamma \) are coefficients that are obtained from experimental data or compiled in chemical databases, and they differ between the VOCs.

Since we know the enthalpy of evaporation and the vapor pressure of a VOC in standard conditions, we may calculate its vapor pressure in the uppermost layer of the stratum corneum. This can be done thanks to the Clausius-Clapeyron equation, which relates the vapor pressure of a substance at two different temperatures:

\[\mathbf{\ln\left(\frac{P_{2}}{P_{1}}\right)} = \] \[\mathbf{- \frac{\Delta H_{\text{vap}}}{R} \left(\frac{1}{T_{2}} - \frac{1}{T_{1}}\right)} \] \[\mathbf{ \Rightarrow e^{\ln\left(\frac{P_{2}}{P_{1}}\right)} = e^{- \frac{\Delta H_{\text{vap}}}{R}\left(\frac{1}{T_{2}} - \frac{1}{T_{1}}\right)} } \] \[\mathbf{\Rightarrow \frac{P_{2}}{P_{1}} = e^{- \frac{\Delta H_{\text{vap}}}{R}\left(\frac{1}{T_{2}} - \frac{1}{T_{1}}\right)} } \] \[ \mathbf{\Rightarrow P_2 = P_1 e^{- \frac{\Delta H_{\text{vap}}}{R}\left(\frac{1}{T_{2}} - \frac{1}{T_{1}}\right)}} \]

2.3 Discussion

• We assume that the Volatile Organic Compounds of interest retain their liquid form in the skin sebum.

• We interpret the existence of the VOCs in the skin sebum as a continuous film. In reality, the VOCs are dissolved in the sebum, and the sebaceous activity differs over the different locations of the skin.

• We assume that the Evaporation rate remains constant during the whole process. In reality, as quantities of the VOC evaporate, and its concentration inside the kit increases, the evaporation rate decreases. However, the diagram that shows this decrease is yet unknown.

• In the literature, there are many formulas that estimate the gas phase mass transfer coefficient. For example, some models include Reynolds Number:

\[ \mathbf{Re = \frac{uL}{\nu}} \]

where L is the evaporation length and $\nu$ is the kinematic viscosity.

This number is related to Nusselt number,

\[ \mathbf{Nu = k_g\frac{L}{D_{a}}} \]

where \(D_{a} \) is the diffusivity of the VOC in air. The equation that relates these numbers comes from Mass Transfer Theory:

\[ \mathbf{Nu = Re^{\gamma } \Rightarrow k_g\frac{L}{D_a} = \left ( \frac{uL}{\nu } \right )^\gamma } \] \[ \mathbf{ \Rightarrow \frac{k_g}{D_a} = \frac{u^\gamma *L^{\gamma -1}}{\nu^\gamma } } \] \[ \mathbf{ \Rightarrow k_g = D_a*L^{\gamma -1}*\left ( \frac{u}{\nu} \right )^\gamma} \]

3. Diffusion

3.1 Evaluating the time required for the kit to remain on the patient's back

In the previous section, we constructed a general formula that estimates the evaporation rate of every VOC of our interest. Plugging every parameter according to the Système International, we end up with a measurement unit for the evaporation rate of \(\text{mol/m}^2\text{s}\)

Since we know the area of the skin that shall be covered by our kit, we can calculate the amount of VOC that evaporates each second:

\[\frac{\Delta n}{\Delta t} = \mathbf{E \cdot A}\]

In order to estimate the amount of time required to fill the volume of the kit with each VOC, we must first determine their desired final concentration. In reality, the concentration itself does not matter for this calculation, as long as the quantity of the VOC is enough for the phage to swell. If this is true, then for every concentration that we decide, the color of the phage that corresponds to a patient with PD will be different from that of a healthy person. Specifically, it is quoted in the literature that the concentration of these VOCs inside the skin sebum of a patient differs from that of a healthy person, and this change is specific for each VOC of interest. Thus, their evaporation rate changes, resulting in different concentrations for the same duration of examination, and in different color outputs. This means that we can begin the experiments with whatever concentration we desire, as long as the phages interact with the VOC. The final result is comparative between the color of the patient and the color of a healthy person, and does not depend on the time we hold the kit on the person's back.

However, since we have to start from some values, we recall the experimental data acquired from Joy Milne, the woman who had congenital hyperosmia and was the inspiration for our project. In reality, the bacteriophage is by far more sensitive as far as the detection of the VOCs is concerned, and this means that the real times are surely lower than the ones calculated.

Mixtures of the candidate compounds (n = 17) were prepared at a range of concentrations (10 \(\mu\)M, 5 \(\mu\)M, 0.5 \(\mu\)M, 0.05 \(\mu\)M, 0.005 \(\mu\)M) and presented to Milne. The results demonstrated she could detect (although not in any systematic order) the whole range of concentrations offered, and a concentration between 0.05 \(\mu\)M and 0.5 \(\mu\)M gave her the best response.

Thus, the desired final concentration of each VOC inside the kit is approximately 0.25 \(\mu\)M.

The volume of our kit is equal to

\[ V = {length} \times {width} \times {height} = \] \[ 5.5 \times 3 \times 2 \times 10^{-6} \, {m}^3 = \] \[33 \times 10^{-6} \, {m}^3 = \] \[ 33 \times 10^{-3} \, {L} \]

We can calculate the amount of VOC required to reach the desired concentration:

\[C = \frac{n}{V} \Rightarrow \] \[ n = C \times V \] \[ \Rightarrow n = 0.25 \, \mu M \times 33 \times 10^{-3} \, L = \] \[ 0.25 \times 10^{-6} \, \frac{\text{mol}}{\text{L}} \times 33 \times 10^{-3} \, \text{L} \] \[\mathbf{\Rightarrow n = 8.25 \times 10^{-9} \, \text{mol}}\]

Taking everything into account, we are now able to calculate the time required to reach the desired amount of every VOC inside the kit.

\[ \begin{Bmatrix} \Delta t_{1} = \frac{n}{E_{1}*A}\\ \\ \Delta t_{2} = \frac{n}{E_{2}*A}\\ \\ \Delta t_{3} = \frac{n}{E_{3}*A}\\ \\ \Delta t_{4} = \frac{n}{E_{4}*A}\\ \end{Bmatrix} \]

Choose

\[ \mathbf{\max\left \{ \Delta t_1, \Delta t_2, \Delta t_3, \Delta t_4 \right \} } \] \[\mathbf{ = \Delta t_{\max}.}\]

The times depend on the 4 VOCs that we examine at a time. By the time we reach \(\Delta t_{max}\), we drag the glass back to its starting position to close the kit, and we remove the latter from the patient's back. We proceed by estimating the time required for the VOCs to evenly spread inside the kit. This means that they will have the desired concentration at the bottom layer of our kit, where the wafers along with the bacteriophages are placed.

Since n and \(A\) are fixed numbers, we may find the correlation between the evaporation rate and \(\Delta t\), which is:

\[ \Delta t = \frac{1}{2} \times 10^{-5} \times \frac{1}{E} \]

Figure 1: Correlation between the fixed Evaporation Rate of a VOC and the required time to fill the kit

This diagram shows the expected value of \(\Delta t\) for a fixed Evaporation Rate, which depends on the VOC. Note that for very small E

\[ \frac{1}{E} \rightarrow +\infty \] \[ \Rightarrow \Delta t \rightarrow +\infty \]

which makes sense, since if the Evaporation rate is extremely slow, the time required for the VOC to fill the kit is extremely large. Vice versa, for large E:

\[ \frac{1}{E} \rightarrow 0^{+} \] \[\Rightarrow \Delta t \rightarrow 0^{+} \]

which also makes sense, since if the Evaporation rate is extremely large, then the required time for the VOC to fill the kit is minimal.

3.2 Evaluating the time required for the VOCs to become evenly spaced inside the kit

By the time we close the kit by dragging the glass, the VOCs will move freely inside the kit, transferring from locations where they have higher concentrations to lower ones.

Figure 2: The physician drags the glass in order to close the kit.

The process continues until the concentration of each VOC is equal in any location of the kit. This equilibrium occurs when the rate of movement due to diffusion is equal in all directions, resulting in a constant concentration profile. We want this process to happen in order to make sure our wafers come into contact with the desired amount of each VOC. In reality, we are calculating the worst-case scenario, as it is almost certain that the wafers will have come into contact with a satisfactory amount of their corresponding VOC, soon before complete diffusion has occurred.

During this phase the screening for PD has been completed and the physician awaits for the results of the kit. The time this process shall take is minimal in proportion to the time required for evaporation. However, if one desires to estimate the former, there are 2 methods:

• Experimentally, after a corresponding statistical analysis: The physician first accumulates a satisfactory amount of results, and then empirically calculates the time required for every VOC to modify the structure of its corresponding bacteriophage. This is the recommended method, as it applies in the real world with no further assumptions.

• By solving the diffusion-convection equation. This method is not recommended for numerous reasons. The diffusion coefficients of the VOCs of our interest remain unrecorded in the literature. Due to the complexity of the equation, there is a vast amount of assumptions that have to be made in order for the former to be solvable. For example, we study the motion of each VOC separately inside the kit, while in reality many compounds exist and change its motion. It is only a predictory formula that can be solved using Computational Fluid Dynamics software. Nonetheless, if one chooses to do so, the diffusion-convection equation is:

\[ \mathbf{\frac{\partial c}{\partial t} = } \] \[ \mathbf{ \nabla \cdot \left ( D\nabla c \right )-\nabla \cdot \left ( uc \right )+R }\]

\[ \mathbf{= D\nabla^2 c - u\nabla\cdot c + R} \]

where D is the diffusion coefficient, u is the airflow velocity, \(\bigtriangledown \cdot\) is the vector operator representing divergence, \( \bigtriangledown ^{2} =\bigtriangledown \cdot \bigtriangledown\) is the Laplacian Operator representing diffusion and R describes possible creation (source) or destruction (sink) of the VOC. Since our kit is closed and the quantity of each voc remains constant,

\[ \mathbf{R = 0} \]

The Laplacian operator is a mathematical operator commonly used in physics and mathematics to describe various physical phenomena, including diffusion. In the n-th dimensional Euclidean Space, and for a twice-differentiable function f, the Laplacian operator is defined as:

\[ \mathbf{\nabla^2 f = \sum_{i=1}^{n} \frac{\partial^2 f}{\partial x_i^2}} \]

over the corresponding coordinates

\[ \mathbf{\left \{ x_1, x_2, . . . ,x_{n-1}, x_n \right \}} \]

As such, the Laplacian operator maps functions of \(C^k\) class to \(C^{k-2}\), being the sum of all the second partial derivatives, and operates on any open set

\[ \mathbf{A \subseteq \mathbb{R}^n_{\text{Euclidean}}} \]

In our case, since we are in a 3 dimensional world,

\[ \mathbf{\bigtriangledown ^2c = \frac{\partial^2c }{\partial x^2} + \frac{\partial^2c }{\partial y^2} + \frac{\partial^2c }{\partial z^2}} \]

Thus, we shall consider the length of our kit as x axis, the width as y, and the height as z. For convenience, we set the height at which our wafers exist as \( \mathbf{max(z)} \), so that the x-y plane which was placed at the back of our patient has a height of \( \mathbf{z=0}\), as displayed in the image.

Figure 3: The 3 axes. Note that z = 0 at the x-y plane that comes into contact with the patient's back.

In order to make the equation easier, one may make the following assumption: by the time we close the kit, all the quantity of each VOC lies on the x-y plane where z = 0, and then starts to spread inside the kit. Of course, this is not true, as the molecules of each VOC will have already spread across the kit, and we will calculate a time bigger than the expected. However, this assumption provides extra information for the equation. Let \(t_{r1}\) be the real time it takes for the VOC '1' to spread evenly inside the kit, and \(t_{a1}\) be the corresponding time, with the assumption mentioned above included. Then we know that

\[ \mathbf{t_{r1} < t_{a1}} \]

The assumption, for \(\mathbf{C(x,y,z,t)}\), can be described as follows:

\[\mathbf{C(x,y,0,0) =C_{max}}\]

and

\[ \mathbf{C(x, y, z, 0) = C_{\text{min}} = 0 } \] \[ \mathbf{ \quad \forall x, y, z, z \not\in \{0\}}\]

where \( \mathbf{C(x,y,z,t)} \) represents the concentration of a VOC at a location with coordinates \( \left \{ x,y,z \right \}\) at a specific time t. We may also add that

\[\mathbf{C(0,y,z,t)=C(x,0,z,t)}\] \[ \mathbf{C(x,y,0,t)= C(x_{max},y,z,t) }\]

\[\mathbf{C(x,y_{max},z,t) = C(x,y,z_{max},t), t\neq 0}\]

to indicate that no quantity of VOCs is being transferred through the walls of our kit.

For a vector field \(\mathbf{F}\), the Divergence operator measures how much the former spreads out or diverges from a given point in space. The divergence of a vector field yields a scalar that represents the rate at which the field flows or diverges away from that point. In our case the divergence operator is defined as

\[\mathbf{\bigtriangledown c \overrightarrow{u} = } \] \[\mathbf{\frac{\partial c \overrightarrow{u}}{\partial x} + \frac{\partial c \overrightarrow{u}}{\partial y}+ \frac{\partial c \overrightarrow{u}}{\partial z}}\]

where \(\bigtriangledown c \overrightarrow{u}\) represents the divergence of the flux of the VOC, and u is the fluid velocity vector.

If the divergence is positive, it indicates that the flux of the VOC is spreading out (diverging) from that point. If it's negative, it suggests that the flux of the VOC is converging toward that point. If the divergence is zero, it implies that there is no net flow of the field at that location.

We finally mention that by the time the molecules of the VOC we examine are evenly spaced inside the kit,

\[ \mathbf{\frac{dC}{dt} = 0 \forall x,y,z}\]

and

\[\mathbf{t_{a1} = } \] \[\mathbf{inf(t: \frac{dC}{dt} = 0 \forall x,y,z)}\]

4. Solving the system of Equations

We remind the reader that we are solving for 4 VOCs simultaneously. Thus our system of equations is:

\[ \begin{Bmatrix} \frac{\partial c_1}{\partial t} = \nabla \cdot \left(D_1 \nabla c_1 \right) - \nabla \cdot \left(uc_1 \right) \\ \\ \frac{\partial c_2}{\partial t} = \nabla \cdot \left(D_2 \nabla c_2 \right) - \nabla \cdot \left(uc_2 \right) \\ \\ \frac{\partial c_3}{\partial t} = \nabla \cdot \left(D_3 \nabla c_3 \right) - \nabla \cdot \left(uc_3 \right) \\ \\ \frac{\partial c_4}{\partial t} = \nabla \cdot \left(D_4 \nabla c_4 \right) - \nabla \cdot \left(uc_4 \right) \end{Bmatrix} \]

Once again, we choose

\[\mathbf{max\left \{ t_{a1},t_{a2},t_{a3},t_{a4} \right \}}\]

and this is the time it takes for the VOCs to evenly spread inside the kit. By adding \(\Delta t_{max}\) to this time, we get the time span from the moment we place the kit on our patient's back, until the final results.

At this point it should be noted that according to literature, when the phages come into contact with their corresponding VOCs, the modification of their structure is almost instantaneous and requires a minimal amount of seconds.

5. Example

We apply the methods mentioned above for octadecanal. We follow this model in order to estimate the required time for the physician to hold the kit on their patient's back.

We know octadecanal's vapor pressure at 25 degrees celsius (\(3.41*10^{-4}\) mm Hg).

We calculate the vapor pressure of octadecanal at 37 degrees Celsius (310 K), using the Clausius-Clapeyron equation.

\[ P_{2\text{oct}} = 0.0009 \, \text{mm Hg} \overset{S.I.}{\rightarrow}0.12 \, \text{Pa} \]

Then, we estimate the Evaporation rate.

The Molecular Weight of octadecanal is \(268.5 g/mol \overset{S.I.}{\rightarrow} 268.5 *10^{-3}kg/mol\). We set u = 0.15 m/s (average value of the range we are testing for) and \(\gamma\) = 0.5 (laminar flow). Substituting, we get:

\[ E_{octadecanal} = 1756 * 10^{-5} *0.123*0.15^{1/2} \] \[ * \frac{\left ( 268.5*10^{-3} \right )^{-1/3}}{8.314*310}\frac{mol}{m^2*s}\]

\[ \approx 5 * 10^{-7} \frac{mol}{m^2s} \]

Thus,

\[ \mathbf{\Delta t_{\text{oct}}= } \] \[\mathbf{ \frac{8.25 \times 10^{-9}}{5 \times 10^{-7} \times 16.5 \times 10^{-4}} \, \text{s}} \]

\[\approx 10 seconds\]

This is the time needed for octadecanal to reach the desired quantity inside the kit. After that, the physician removes the kit from the patient's back, closes the former by dragging the glass back to its starting position, and then diffusion happens. Since diffusion is occurring in atmospheric air inside a kit that has a small volume, we assume that the process should not take more than a few seconds. After diffusion has occurred, we will be sure that octadecanal's wafer has come into contact with the desired VOC. This means that the engineered phages that remain encrusted on the wafer will swell, changing their color. This structural modification is almost instantaneous, according to the literature.

Following the same steps, we have determined the corresponding time for other VOCs of interest.

\[ \begin{align*} E_{\text{eicosane}} & : 1.2 \times 10^{-8} \, \frac{\text{mol}}{m^2s} \\ \Delta t_{\text{eicosane}} & \approx 416 \, \text{seconds} \\ \\ E_{\text{perillaldehyde}} & : 6.02 \times 10^{-5} \, \frac{\text{mol}}{m^2s} \\ \Delta t_{\text{perillaldehyde}} & \approx 0.08 \, \text{seconds} \\ \\ E_{\text{hippuric}} & : 1.545 \times 10^{-9} \, \frac{\text{mol}}{m^2s} \\ \Delta t_{\text{hippuric}} & = 3230 \, \text{seconds} \\ \end{align*} \]

We finally remind the reader that the concentration used to obtain these results was taken from the experimental data by Joy Milne. The M13 Bacteriophage is extremely sensitive as far as detecting the desired VOCs is concerned, thus we can obtain the change of color that we want with even lower concentrations of VOCs. This will vastly reduce the wait time for VOCs like hippuric acid. The concentration is used in this model only to calculate the estimated time of examination. The final result (color change) is comparative between a healthy and a PD person, and it will be different independently of the selected concentration.