Team:TU Delft/Model/Q5

iGEM TU Delft

Modeling

Question 5:

How does the polysilicate layer covered cell focus the light?

COMSOL modeling

Short intro about COMSOL Multiphysics:

The software COMSOL Multiphysics was used to model the electromagnetic field’s interaction with our structure. COMSOL Multiphysics is a CAE software package that can be used to model any physics based system and the interaction of different physics (COMSOL Multiphysics®, 2016). For this project the RF Module was used; this module is the best fit for the dimensions of our structure (micrometers) and the intended wavelength (in the visible spectrum). The figure below demonstrates the different options for different structure size and wavelengths (COMSOL, 2013).

COMSOL modules
Figure 1: Best COMSOL modules to use relevant to the Object size and wavelength (COMSOL, 2013).

3D model – Sphere

For the Electromagnetic field–structure interaction there are two models built. The first and simplest assumes that the shape of our structure is a simple sphere. Using the simple sphere model we have symmetry allowing to model only a quarter of the sphere and we need smaller domain compared to a rod shaped structure resulting in a computational less expensive model and faster conversions. Also the case of spherical lenses is better design due to their independence from their orientation to the incoming radiation compared to the bacterial shaped lenses.

Furthermore parametric method of modelling was used, meaning that the most important parameters for the model were defined and then used for the model design. The parameters used for this model is the radius of the sphere, the wavelength, wavenumber and frequency of the intended light, the thickness of the silicate layer, air layer and Perfectly Matched Layer, the intensity of the intended electromagnetic field, and the material parameters epsilon defined later. The values of those parameters are shown in table 1 below.

Table 1: Values of parameters used in the model

ParameterValueDescription
r05*10-7 [m]Radius of the cell
lda5*107 [m]Wavelength
k01.2566*107 [1/m]Wavenumber in vacuum
f05.9958*1014 [1/s]Frequency
t_medium2.5*10-7 [m]Thickness of air layer
t_pml2.5*10-7 [m]Thickness of Perfectly Matched Layer
t_sil8*10-8 [m]Thickness of silicate layer
E01 [V/m]Intended electromagnetic field

According to the aforementioned parameters a three dimensional model of a layered sphere, representing our structure, was created in COMSOL Multiphysics®. The 3D model can be seen in Figure 2. In that model the inner part is the r0, representing the cell, the first layer is the silicate layer covering the cell, then the other two layers are the medium with the outermost representing the world further away from the structure. It is important to note that the 3D model shown in Figure 2 is only a quarter of the actual structure but due to symmetry it is possible to use only that part to decrease the computational cost.

3D sphere in small domain
Figure 2: 3D design of the sphere in a small domain.

Materials

Next the materials need to be determined. The RF module uses three important material parameters for its calculations: relative permeability ( \( \mu_r \) ), electrical conductivity ( \( \sigma \) ) and relative permittivity (\( \epsilon_r \)). The relative permeability is equal or almost equal to unity for most real materials for the optical frequency range that concern us (visible spectrum of the EM field) (Mcintyre and Aspnes, 1971). The values for the electrical conductivity were obtained from the material library of COMSOL Multiphysics® and from literature and the relative permittivity can be calculated from the refractive index (n) using the following formulas(Griffiths, 1999):

$$\epsilon = \epsilon' – j \epsilon'' $$ $$ \epsilon' = \frac{n^2 – k^2}{\mu} = n^2 – k^2$$ $$\epsilon'' = 2 \times n \times \frac{k}{\mu} = 2 \times n \times k$$

So it is:

$$\epsilon = n^2 – k^2 – j \times 2 \times n \times k $$

Here n is the Real part of the refractive index and k is the Imaginary part. Because we assume that we have non absorbing materials and thus the complex part of the refractive index is zero ( \( k=0 \) ) the relative permittivity can be calculated from the refractive index as:

$$\epsilon = n^2$$

Concluding the material parameters used for this model are summarized in Table 2 below. The refractive index of water is 1.33 (Daimon & Masumura, 2007) and of the cell 1.401 (Jericho, Kreuzer, Kanka, & Riesenberg, 2012) and the relative permittivity of both is calculated using the aforementioned formula. The same method was used to calculate the relative permittivity of Tin dioxide, with refractive index between 2.33 and 2.8 for 550 nm (Baco, Chik, & Md. Yassin, 2012) the relative permittivity is 6.58.

Table 2: Material parameters used in the model

Parameters/MaterialsMedium (water)Glass LayerTin dioxide (SnO2)Cell
Relative permeability (μr)1111
Electrical conductivity (σ)0.05 [S/m]10-14 [S/m]0.025 (Banyamin, et al., 2014)0.48 [S/m] (Castellarnau, et al., 2006)
Relative permittivity (εr)1.772.096.581.96

Electromagnetic wave

The intended electric field is \(E0*\exp(-j*k0*x)\) and passes through the whole structure. Figure 3 below demonstrates how EM radiation propagates in space, we have set the propagation direction as x, the electric field oscillation as z and of course the magnetic field oscillation as y.

pherical model
Figure 3: Electromagnetic radiation propagation.

Mesh

A very important part of finite element modelling (FEM) is the meshing of the design. The important part is that there should be enough nodes that the structure is well represented from the FEM but not so many that the system runs out of memory and the simulation never finishes. The meshing is very important for our model because we have a very thin silicate layer between the cell and the air. We cannot model it with coarse element size as there wouldn’t be enough detail and the simulation will never converge and on the other hand we can’t mesh the whole structure finely because it is unnecessary. Keep in mind that when changing the mesh, the number of the Degrees of Freedom of the model increase in the number of 3 because we are using a 3 dimensional domain. The general rule of thumb for meshing in RF simulations is to use as maximum element size about one tenth of the wavelength, in our case this is about 50 nm. The meshed structure can be seen in Figure 3.

Meshed spherical model
Figure 4: Meshed Structure.

3D - Rod Shaped model

Additionally to the spherical model, a rod shaped model that resembles the shape of bacteria closer was used. The modeling method, parameters and materials used to create the rod shaped structure are the same as the spherical one. The way this was modeled is with a layered cylinder and two half spheres in each end of the cylinder. The 3D model of the rod can be seen in the figure below. The length of the middle part was set to \({0.5} \mu m\)

Rod Model
Figure 5: 3D design of the rod shaped structure.

The meshing was created with the same rule as well, maximum mesh size selected was again 50 nm for this structure.

Meshed rod model
Figure 6: Meshed rod shaped structure.

Studies performed

First simulations with small domain

The aforementioned models were used to predict the behavior of light when it meets our biolenses. The first simulation was that of a circular cell of diameter \( 0.5 \mu m \) covered with a thin film of \( 80 nm\) silica. The domain selected was water and the domain was circular for better use of symmetry. The domain in this study was very small at \( 1 \mu m\). The reason a small domain was selected is the small computational time during the troubleshooting period. The first results of those simulations can be seen below.

”Z
Figure 7: Electric field, z component.

Figure 7 shows the z component of the electric field. It can be seen that some focusing is present in the structure. This can be seen better in figure 8, the normalized electric field. The focusing here can be seen in the edge of the cell. One reason that we see that can be that the focal point of the lens is further than the \(1 \mu m\) of this simulation and it is somewhere in the Perfectly Matched Layer (PML) where we can’t see it. In order to investigate if this is the reason we can’t see a focal point additional simulations were performed with larger domain. Here needs to be noted again that while the domain is enlarged the number of nodes increase in the power of 3 so the increase of domain for example to \(1 mm\) is impossible at least with the equipment we have.

One additional comment on this model is that the intended electric field can be seen propagating with no change except from the part that interacts with the sphere and a little around it due to some scattering.

Normalized electric field.
Figure 8: Electric field normalized.

Silica covered cells

Since the size of the PHB granules that can be grown inside an E. coli bacterium was expected to be smaller than the required size for lasing, we investigated the option of making a biolaser out of the entire cell. Therefore we encapsulated the cell with a layer of biosilica and tin dioxide. The gain medium in this method is provided by fluorophores which we expressed in the cytosol of the cell. To determine the minimal size in this case the calculation was a bit more tricky than for the PGB granules, since we had to take the layer into account and the optical path has a star-like shape (figure 3A) .

Figure 7: The path of light in a cell covered with a layer of biosilica or tin dioxide. The optical path has a star like shape.

The critical angle can again be computed by equation 3, where \(n_1\) is the refractive index of the layer \(n_f\), and \(n_2\) is the refractive index of the buffer outside the cell \(n_b\). When the angle is slightly larger than the critical angle we have total internal reflection where the outgoing angle is equal to the incoming angle, therefore we may set \(\alpha = \theta_c\). In this model we will neglect the curvature of the surface which means that we can also set \(\theta_i=\alpha\). Using Snell’s Law (equation 2) again we can compute the outgoing angle \(\theta_0\) (equation 10) where \(n_f\) is the refractive index of the layer and \(n_c\) the refractive index of the cytosol.

$$\theta_0 = asin\Big(\frac{n_f}{n_c} sin(\alpha)\Big) $$

From \(\theta_0\) we can easily determine the angle \( \phi \) since all the angles in a triangle add up to \(\pi\) (equation 11).

$$\phi = \pi-2 \theta_0 = \pi-2\arcsin\Big(\frac{n_f}{n_c} sin(\alpha)\Big)$$

Here the path of the light is star-shaped and the \( OPL=k_{sides}\lambda\). Furthermore we will assume that the total path length per side is minimally \(\lambda\) and that one side is taken as the path from point A to D. The optical path length per side is given as \(OPL_{side}=2\cdot x+l\geq\lambda\). \(X\) and \(l\) can be determined as in equation 12. From equation 12 the minimal radius can be determined for the criteria where \(OPL_{sides} = \lambda\) (equation 13).

$$\left.\begin{matrix} OPL_{sides} =2x+l\geq \lambda \\ x=\frac{t}{\cos\alpha}\\ l =2R\cos\theta_o \end{matrix}\right\} OPL_{sides} = \lambda = \frac{2t}{\cos\alpha}+2R \cos\theta_o =\frac{2t}{\cos\alpha}+2R \sqrt{1-(\frac{n_f}{n_c}\sin\alpha)^2}\geq\lambda $$ $$ R_{min} = \frac{\lambda-\frac{2t}{\cos\alpha}}{2\cos\theta_o}$$

Using the minimal radius of R we can determine the angle \(\psi\) as in equation 14.

$$ \psi=2\arcsin(\frac{t \tan\alpha}{R_{min}})$$

Using \(\psi\) and \(\phi\) we can determine the minimal number of sides of the optical path (equation 15).

$$k_{sides} = \frac{2\pi}{\psi+\phi}$$

In equation 15, \(k_{sides}\) has to be an integer to have a closed optical path. Taken together the constraints as in equation 16 we get an equation for \(\phi\) and \(\psi\) as a function of \(\alpha\).

$$\left.\begin{matrix} \frac{2\pi}{k_{sides}}=\phi+\psi\\ \phi = \pi-2\arcsin(\frac{n_f}{n_c}\sin\alpha)\\ \psi = 2\arcsin(\frac{t\tan\alpha}{R_{min}})\\ OPL_{sides}=2x+l \rightarrow R_{min} = \frac{\lambda-\frac{2t}{\cos\alpha}}{2\sqrt{1-(\frac{n_f}{n_c}\sin\alpha)^2}} \end{matrix}\right\} \frac{2\pi}{k_{sides}}=\pi-2\arcsin(\frac{n_f}{n_c}\sin\alpha)+2\arcsin(\frac{2t\sin\alpha\sqrt{1-(\frac{n_f}{n_c}\sin\alpha)^2}}{\lambda\cos\alpha-2t})$$

Solving equation 16 for \(\alpha\) and using this in the equation for the radius results in the minimal radius required for WGM in a cell covered with a layer of biosilica or tin dioxide.

Results

Using the model described above we determined the minimal radius of a cell covered with polysilica and tin dioxide to create a cavity that allows for whispering gallery modes. The refractive index of polysilica is 1.47 (Liang et al., 2007). The refractive index of tin dioxide is 2 (Pradyot, 2003). For both cases we calculated the minimal radius assuming the thickness of the layer to be 50 nm and a wavelength of 509 nm, the emission wavelength of GFP. From our model we found a minimal diameter of \(1.2 \mu m\) and \(1.6 \mu m\) (figure 8) for polysilica and tin dioxide respectively. These sizes are comparable to the size of a bacterium, so whispering gallery modes should in principle be possible inside a covered cell. However, we made a number of simplifications in our model that will in practice probably prevent the system from lasing at this minimum scale. First, we did not take the polarization of the light into account. Second, the length scales we work with here are only about twice the wavelength, and the thickness of the boundary layer is much smaller than the wavelength. A more accurate description of what is happening within the cell would require us to work in the thin layer limit (known as Mie theory) as we do in Q4. There we calculate the quality factor for the cavity, which indicates how well the light can be trapped within the structure. Third, in this model we only determined whether light waves will fit into the cavity, in Q3 we will also determine whether lasing can take place with the size of the cell and the reflectivity of polysilica and tin dioxide.

Figure 8: (A) For a cell covered in a 50 nm polysilica the minimal diameter is \(1.2\mu m\). (B) For a cell covered in tin dioxide the minimal diameter is \(1.6\mu m\).
  1. Aizenberg, J., Sundar, V. C., Yablon, A. D., Weaver, J. C., & Chen, G. (2004). Biological glass fibers: Correlation between optical and structural properties. Proceedings of the National Academy of Sciences of the United States of America, 101(10), 3358–3363. https://doi.org/DOI 10.1073/pnas.0307843101
  2. Hecht, E. (2001). Optics 4th edition. Optics 4th Edition by Eugene Hecht Reading MA Addison Wesley Publishing Company 2001. https://doi.org/10.1119/1.3274347
  3. Huglin, M. B., & Radwan, M. A. (1991). Behaviour of poly(β-hydroxybutyric acid) in dilute solution. Polymer International, 24(2), 119–123. https://doi.org/10.1002/pi.4990240210
  4. Humar, M., & Yun, S. H. (2015). Intracellular microlasers. Nature Photonics, 9(9), 572–576. https://doi.org/10.1038/nphoton.2015.129
  5. Kitamura, S., & Doi, Y. (1994). Staining method of poly(3-hydroxyalkanoic acids) producing bacteria by nile blue. Biotechnology Techniques, 8(5), 345–350. https://doi.org/10.1007/BF02428979
  6. Liang, X. J., Liu, A. Q., Lim, C. S., Ayi, T. C., & Yap, P. H. (2007). Determining refractive index of single living cell using an integrated microchip. Sensors and Actuators A: Physical, 133(2), 349-354.
  7. Patnaik, P. (2003). Handbook of Inorganic Chemicals. Ebook. Retrieved from ftp://pvictor.homeftp.net/public/Sci_Library/Chem Library/Handbooks/Patnaik P. Handbook of inorganic chemicals (MGH, 2003)(T)(1125s).pdf
  8. Strutt, J. W., & Rayleigh, B. (1877). The Theory of Sound. Nature (Vol. 1). https://doi.org/10.1038/058121a0
  9. Wyatt, P. J. (1968). Differential Light Scattering: a Physical Method for Identifying Living Bacterial Cells. Applied Optics, 7(10), 1879. https://doi.org/10.1364/AO.7.001879