SOCR ≫ BPAD Website ≫ BPAD GitHub ≫

Earlier, in the mathematical foundations chapter we discussed the concepts of scalars, vectors, matrices, and tensors, as well as some of the corresponding (tensor) operations, e.g., scalar multiplication, vector product, and matrix multiplication.

1 Distances and Sizes

In its most general form, the term distance is defined as the amount of “space” separating two objects, states, or even abstract concepts, where space is the union of all possible objects or states. In geometric settings, distance between two points represents a numerical measurement reflecting the size of space separating the points. In physics, a distance measure characterizes the physical length between particles or objects. For a pair of objects, think about two points \(A\) and \(B\), the distance between them may be denoted by \(dist(A,B)= |AB|\equiv |BA|=dist(B,A)\). The mathematical description of a distance function (metric) describes proximity of elements of some given space - nearby and far-apart objects correspond to small and large distances, respectively.

Examples of distance metrics, i.e., space-metric pairs, include:

  • The Euclidean distance between two points in the 2D plane. Pick any pair of points \(P_1=(x_1, y_1),\ P_2=(x_2, y_2)\in \mathbb{R}^2\equiv \mathbb{C}\), then, the distance between them is

\[d(P_1,P_2)=\sqrt {(\Delta x)^2+(\Delta y)^2} =\sqrt {(x_2-x_1)^2+(y_2-y_1)^2}.\] This formulation of a distance in \(\mathbb{R}^2\) is derived by applying the Pythagorean theorem. The result is a 2D flat (Euclidean) metric-space \(\left(\mathbb{R}^2, d\right)\). Note that the distance \(d\) leads to the arc-length distance quantifying the distance between pairs of points along a curve embedded in \(\mathbb{R}^2\).

There are many other distance metrics that can be defined to construct \(n\)-dimensional Euclidean and non-Euclidean metric-spaces. For instance the Euclidean metric-space \(\mathbb{R}^n\) can be equipped with some of the following distance measures \(d_p(X,Y):\mathbb{R}^n\times \mathbb{R}^n \to \mathbb{R},\ \forall p\in \mathbb{R}_{\ge 1}\) (Minkowski distance of order \(p\), or simply \(p\)-norm):

\[d_p(X,Y) =\left(\sum_{i=1}^{n}\left|x_i-y_i\right|^{p}\right)^{\frac{1}{p}}\equiv ||X-Y||_p.\] However, the infinity norm distance, or Chebyshev distance, provides yet another way to track object separation in \(\mathbb{R}^n\): \[d_{\infty}(X,Y)=\lim _{p\to \infty}\left(\sum_{i=1}^{n}\left|x_{i}-y_{i}\right|^{p}\right)^{\frac{1}{p}} \equiv \max \left(|x_1-y_1|,|x_2-y_2|,\cdots ,|x_n-y_n|\right)\equiv ||X-Y||_{\infty}.\]

In biophysical problems, the experimental space, e.g., \(\mathbb{R}^n\), is often equipped with a Euclidean distance of the corresponding dimension, e.g., \(d_{2,\mathbb{R}^n}(\cdot,\cdot)\). This choice provides a natural metric-space to study the phenomenon where length/distances are invariant with respect to rigid-body transformations, e.g., rotation and offset/shift.

All distance measures must satisfy three metric conditions. On the space \(M\), a distance function \(d\) is a map \(d(\cdot,\cdot):M\times M \to \mathbb{R}^+\) must obey the following conditions:

  • \(d(x,y) \geq 0,\ \forall x,y\in M\), with \(d(x,y) = 0 \iff x = y\)
  • \(d(x,y) = d(y,x),\ \forall x,y\in M\)
  • \(d(x,z) \leq d(x,y) + d(y,z),\ \forall x,y,z\in M\).

In the simple situation where \(M\equiv \mathbb{R}\) is the space of the real scalars, the natural definition of the distance between a pair real numbers \(x, y\in\mathbb{R}\) is \(d(x,y) = |x − y|\), the absolute-value of the numerical difference between the numbers. Check that this distance satisfies the three metric conditions. An alternative distance formulation on the same state-space is \(\tilde{d}(x,y) = \begin{cases} 0 & x=y \\ 1 & x\not= y \end{cases}\), which defines a different metric corresponding to another topological configuration of the real line called discrete topology where all reals are dispersed and can never be arbitrarily close to each other (say within \(0.1\) units from one another).

The size of a space, e.g., state space of a biophysics experiment, is measured using a properly defined metric on the space. For instance, length (1D space), area (2D space), volume (3D space), and and hyper-volume (\(n\)D space) represents commonly used quantitative measures to capture the size of a neighborhood, or a subspace/subset, in the original space. Length naturally corresponds to sizes of 1D linear segments, whereas area and volume capture the size of a closed region in a 2D region in the plane and a 3D solid, respectively. For instance, using the classical Euclidean distance, the area of a square in 2D, with side of \(1\ \text{unit}\), is equal to \(1\times 1=1\ \text{unit}^2\) and the volume of a ball of radius \(r\) anywhere in 3D is \(\frac{4}{3}\ \pi\ r^3\).

The FSU Galactic-to-Microscopic animation shows the Universal scale from a distant cosmic view of the Milky Way Galaxy, from a distant viewpoint 10 million light years away from Earth (\(10^{23}\) meters), down to the peculiar world of the quark particles (\(10^{-16}\) meters) by orders of magnitude magnification (\(\times10\)). This universal view starts with relativistic distances and goes through the Newtonian and human scales into the enigmatic microscopic world of biological cells, DNA, electrons, and protons.

Biophysical systems study objects spanning a wide range of sizes, scales, shapes, and forms. In the metric system, the meter is the universal unit of length. Life diversity increases with lowering length scale, i.e., the variety and dynamics of microorganisms at lower scales are much more pronounced than their complex-life-form counterparts at larger scales.

Human vision allows us to observe experiments with naked eyes for scales from around \(1\ km\equiv 10^3\ m\) to about \(1\ mm \equiv 10^{−3}\ m\). Of course, we use instruments to enhance our senses, peer outside of this limited range. Advanced hardware, software, and artificial intelligence allow us to study detailed structure of distant galaxies light-years away (\(\ge 10^{16}\ m\)) or tiny objects \(1\ \mu m\equiv 10^{−6}\ m\) (micrometer, micron), \(10^{−10}\ m\) (nano=-meter), etc.

2 Models and Process Modeling

There are a number of complementary formulations of models based on the scientific domain, the application construct, and the representation framework. Some examples of commonly used models include the following:

  • Generic model: a simplified representation of an object, a complex system, or a concept. Think of a maquette of something like a building, a structure, a city, or even a cosmic system of celestial bodies.

  • Physical model: a minimal, idealized, and real representation of an object, e.g., a toy car, that illustrates its key characteristics and properties, e.g., car parts, motion dynamics, functionality, limitations, etc.

  • Abstract mathematical model: a conceptual mathematical expression, typically derived as a set of equations, describing patterns, structures, shapes, forms, relations, inter-dependencies, which permit symbolic, arithmetic, numeric, or functional operations.

  • Design and art model: a template, product, mold, or object that is itself incomplete but shows the important characteristics of the corresponding real object it represents.

  • Biophysical model: a model organism, cell-line, species, or a testable maquette that can be used for in vitro studies aiming to understand a given biological system, medical phenomena, or a mechanistic biophysical principle.

There are a number of other types of models, including conceptual, algorithmic, econometric, statistical, logical, software, axiomatic, performance and acting, societal, and political models.

Let’s consider a couple of relevant biophysical models.

2.1 Harmonic oscillator model

Suppose we try to understand the dynamical mechanism behind the motion of an object of mass \(m\) attached to one end of a spring whose other end is held fixed. Under the tug between gravity and the spring forces, when the object is disturbed from its steady-state (at rest equilibrium), the motion of the object is expected to oscillate along a straight line about its equilibrium position.

Suppose the oscillations take place in the direction of the vertical \(z\)-axis and the position of the origin of the coordinate system is fixed at the equilibrium position. According to the Hooke’s law of physics, the spring exerts a force \(F\) on the object of mass \(m\), which drives the object towards its equilibrium position. Under certain conditions, the magnitude of the force is proportional to the distance from the object to its equilibrium position. To illustrate a mathematical model of the physical system, assume we are given a spring with a corresponding spring-constant \(\kappa \in\mathbb{R}^+\), which depends on the spring material, construction protocol, and tightness. Then the Hooke’s law is

\[ F = −\kappa x\vec{i},\]

where \(F\) is the force acting on the object, \(\vec{i}\) is a unitary vector in the positive direction of the motion, e.g., vertical, \(\kappa\) is the spring constant, and \(x\) is the displacement amount from the point of equilibrium.

The Newton’s second law yields that the total force \(F\) acting on a mass \(m\) is proportional to the acceleration \(a=x''\) of the object resulting from the force, and the mass of the object represents this constant of proportionality:

\[F = m\ a.\]

For simplicity of the model formulation and ease of its interpretation, we can assume that there are no other forces affecting the spring dynamics, i.e., no other actions on the mass, no environmental frictions, no mass-loss, no gravity, no confounding wind or electromagnetic fields.

Then, at each time instant \(t\) (the dynamics parameter), equating the right-hand-sized of the Hook’s and Newton’s laws yields the following ordinary differential equation model describing the spring motion dynamics

\[m\ x''(t) = −\kappa x(t),\]

where \(m\) is the mass of the object and \(x''(t)\) is the acceleration, i.e., the second-order time-derivative of the displacement \(x(t)\), which represents the object spatial position at time \(t\). This model can also be expressed as

\[x''(t) +\omega^2 x(t)=0,\]

where \(\omega = \sqrt{\frac{\kappa}{m}}\).

This equation is called the harmonic oscillator equation model. It is a simple homogeneous, second order, linear differential equation with constant coefficients whose general (ansatz) solution expressing the spatial displacement at a given time point is:

\[x(t) = A \cos (\omega t + \varphi), \] where \(A\) is the oscillation amplitude, \(\omega\) is the natural frequency, and \(\varphi\) is the phase of the motion. Validate that for \(A=\frac{\kappa}{\omega^2}\), \(m\ x''(t)\equiv −\kappa x(t)\): \[m\ x''(t) = \omega (-A\sin(\omega t + \varphi))'= -\omega^2 A \cos (\omega t + \varphi) \equiv −\kappa x(t).\]

Of course, we can relax some of the strong assumptions we made above to account for effects of other forces (damping, friction, or driving forces).

2.2 Color vision model

Vision is one of the main five senses used to perceive an ambient environment. In addition to tracking shades, shapes and dynamical changes, vision can be significantly augmented by an ability to detect existence of color. Color is associated with the corresponding wavelengths of the object light reflections in the visible electromagnetic spectrum. Eyes receive light wavelengths as colors, such as hues of pure or mixtures of red, orange, yellow, green, blue, and violet, spanning the rainbow of colors decomposing white light into different wavelength bands. The retina in the human eye has two types of photoreceptors (light-sensing cells) - rods and cones.

The eye has about \(120M\) rods that that are located peripherally and are \(\sim 1,000\) times more sensitive than cones. Rods are responsible for peripheral vision, low-light vision, and motion detection, but less sensitive to color. the \(\sim 6M\) cones are most concentrated in the central region of the retina (fovea) and are responsible for high resolution vision. The three types of cones are sensitive to different ranges of wavelengths corresponding to red, green, and blue spectra (RGB spectra). A wide range of hues are perceived by mixing appropriate ratios of inputs form the base three colors, e.g., yellow is a blend of red and green, and white is a mixture of all three base colors. The figure below shows a Venn diagram mixing red, green and blue colors.

# install.packages("VennDiagram")
library(plotly)
library(VennDiagram)

# Define pair/triple set intersections
n12 = 320+70; n23 = 190+70;  n13 = 60+70; n123 = 70

# Generate Venn diagram but don;t plot it yet # ind = FALSE
plt <- draw.triple.venn(area1 = 400+60+70+320, area2 = 800+190+320+70, 
                        area3 = 3200+190+60+70, n12 = 320+70, n23 = 190+70, 
                        n13 = 60+70, n123 = 70, col="gray", 
                        category = c("Red", "Green", "Blue"), 
                        lty = 1, lwd=4, fill = c("Red", "Green", "blue4"), 
                        alpha = c(0.5, 0.3, 0.4), label.col=c("Red",
                          "Yellow","Green", "Purple", "White", "Brown", "Blue"),
                        cex=2, ind = FALSE)

# this plot object (plt) is complex - a list of 8 elements, each element
# contains a polygon, e.g., # i=2; plot(x=plt[[i]]$x, y=plt[[i]]$y)

text <- data.frame(x=c(0.2, 0.5, 0.8, 0.3, 0.5, 0.7, 0.5), 
                   y=c(0.7, 0.8, 0.7, 0.4, 0.5, 0.4, 0.1), 
                   value=c(400, 320, 800, 60, 70, 190, 3200))

plot_ly(type="scatter", mode="lines", x=plt[[1]]$x, y=plt[[1]]$y, 
        name="Set 1", opacity=0.5, line=list(color = 'red', width = 4),
        fill = 'toself', fillcolor = 'pink', hoveron='points') %>%
  add_trace(x=plt[[2]]$x, y=plt[[2]]$y, name="Set 2", opacity=0.5, 
        line=list(color = 'green', width = 4), fill='toself', 
        fillcolor='green', hoveron='points') %>%
  add_trace(x=plt[[3]]$x, y=plt[[3]]$y, name="Set 3", opacity=0.5, 
        line=list(color = 'blue', width = 4), fill='toself', 
        fillcolor='blue', hoveron='points') %>%
  # Extra sets are added as additional traces
  # add_trace(x=plt[[4]]$x, y=plt[[4]]$y, name="4") %>%
  # add_trace(x=plt[[5]]$x, y=plt[[5]]$y, name="5") %>%
  # add_trace(x=plt[[6]]$x, y=plt[[6]]$y, name="6") %>%
  # add regional intersection value labels, as annotations
  layout(legend = list(orientation = 'h'), title="Dynamic Venn Diagram",
  annotations=list(x=text$x, y=text$y, text=text$value))

Try to save the following R code as a stand-alone RShiny app (RGB_ColorApp.R) and run it to see how RGB mixtures of primary colors constitute a wide spectrum of visible colors.

library(plotly)
library(VennDiagram)
library(shiny)
library(shinyjs)
library(shinydashboard)
library(plothelper)
library(colourpicker)

shinyApp(
  ui = fluidPage(
    fluidRow(
      column(3, colourpicker::colourInput("color1", "Set 1 Color", "red")),
      column(3, colourpicker::colourInput("color2", "Set 2 Color", "green")),
      column(3, colourpicker::colourInput("color3", "Set 3 Color", "blue")),
      column(6,plotlyOutput("VennDiagram"))
    )
  ),
  
  server = function(input, output) {
    output$VennDiagram <- renderPlotly({
      
      plot_ly(type="scatter", mode="lines", x=plt[[1]]$x, y=plt[[1]]$y, 
              name="Set 1", opacity=0.5, line=list(color = input$color1, width = 4),
              fill = 'toself', fillcolor = input$color1, hoveron='points') %>%
        add_trace(x=plt[[2]]$x, y=plt[[2]]$y, name="Set 2", opacity=0.5, 
                  line=list(color = input$color2, width = 4), fill='toself', 
                  fillcolor=input$color2, hoveron='points') %>%
        add_trace(x=plt[[3]]$x, y=plt[[3]]$y, name="Set 3", opacity=0.5, 
                  line=list(color = input$color3, width = 4), fill='toself', 
                  fillcolor=input$color3, hoveron='points') %>%
        layout(legend = list(orientation = 'h'), title="Dynamic RGB Color Venn Diagram",
               annotations=list(x=text$x, y=text$y, text=text$value))
      # insert code
    })
  }
)

In practice, this model of color vision represents a simple foundation for interpreting a much more complex reality of color vision where various different color combinations can produce the same hue. For instance, yellow can be sensed with a combination of red and green or as white light with removed violet.

The figure below shows the color perception of color objects illuminated with white light (complete wavelength spectrum). In this case, a blue object absorbs all wavelengths except the in the blue spectrum.

library(plotly)
library(plotly)
## Light rays (sources)
# Blue
blueEdges <- list(x = c(0.2, 1), y = c(0.2, 1), z = c(0.8, 0), name="Blue", showscale=F,
  line = list(color = "lightblue",  width = 15), mode = "lines", opacity = 1, type = "scatter3d")
blueEdges1 <- list(x = c(1, 1.5), y = c(1,1.5), z = c(0, 0.5), name="Blue1", showscale=F,
  line = list(color = "lightblue",  width = 10), mode = "lines", opacity = 1, type = "scatter3d")
# Green
greenEdges <- list(x = c(0.4, 1.2), y = c(0.0, 0.8), z = c(0.8, 0), name="Green", showscale=F,
  line = list(color = "green",  width = 15), mode = "lines", opacity = 1, type = "scatter3d")
greenEdges1 <- list(x = c(1.2, 1.7), y = c(0.8, 1.3), z = c(0, 0.5), name="Green1", showscale=F,
  line = list(color = "green",  width = 1, opacity = 0.1), mode = "lines", opacity = 0.1, type = "scatter3d")
# Red
redEdges <- list(x = c(0.0, 0.8), y = c(0.4, 1.2), z = c(0.8, 0), name="Red", showscale=F,
  line = list(color = "red",  width = 15), mode = "lines", opacity = 1, type = "scatter3d")
redEdges1 <- list(x = c(0.8, 1.3), y = c(1.2, 1.7), z = c(0, 0.5), name="Red1", showscale=F,
  line = list(color = rgb(0.5,0,0,0.5),  width = 1), mode = "lines", type = "scatter3d")

# Scene color restyle buttons
updateColorMenus = list(list(name = "Reflection<br>Surface Color"), 
                   list(type = "buttons", y = 0.8,
                        buttons = list(list(method = "restyle",
                                         args=list(list(facecolor=list(c("red","red")),
                                                 visible=c(T,T,F,T,T,T,T,F,F,T,F,T,T))),
                                            label = "red"),
                                       list(method = "restyle",
                                         args=list(list(facecolor=list(c("green","green")),
                                                 visible=c(T,T,F,T,T,T,T,F,T,F,T,F,T))),
                                            label = "green"),
                                       list(method = "restyle",
                                        args=list(list(facecolor=list(c("blue","blue")),
                                                 visible=c(T,T,T,T,T,T,T,T,F,F,F,F,F))),
                                            label = "blue"),
                                       list(method = "restyle",
                                        args=list(list(facecolor=list(c("black","black")),
                                                 visible=c(T,T,F,T,T,T,T,F,F,F,F,F,T))),
                                            label = "black")
                                       )))

# plot_ly(name="Blue Object", x=c(0,2,2,0), y=c(0,0,2,2), z=c(0,0,0,0), # Points
#         i = c(1, 0), j = c(2,1), k = c(3,3), # faces
#         type = 'mesh3d', showscale=F, facecolor = c("red","red")) %>%
#   layout(legend = list(orientation = 'h'), 
#          title="Blue Objects Absorb all Wavelents and Reflect Blue",
#          scene = list(aspectratio = list(x = 1,y = 1,z = 1)),
#          updatemenus = updatemenus)

## Plot 3D Scene
plot_ly(name="Object Surface", x=c(0,2,2,0), y=c(0,0,2,2), z=c(0,0,0,0), # Points
        i = c(1, 0), j = c(2,1), k = c(3,3), # faces
        type = 'mesh3d', showscale=F, facecolor=c("blue","blue")) %>% 
  add_trace(x=blueEdges$x, y=blueEdges$y, z=blueEdges$z, type="scatter3d", 
            mode="lines", name="BlueLight", visible=FALSE,
            line = list(color=blueEdges$line$color, width=blueEdges$line$width), showscale=F) %>%
  add_trace(x=blueEdges1$x, y=blueEdges1$y, z=blueEdges1$z, type="scatter3d", mode="lines", name="BlueReflect",
            line = list(color=blueEdges1$line$color, width=blueEdges1$line$width), showscale=F) %>%
  add_trace(x=greenEdges$x, y=greenEdges$y, z=greenEdges$z, type="scatter3d", mode="lines", name="GreenLight",
            line = list(color=greenEdges$line$color, width=greenEdges$line$width), showscale=F) %>%
  add_trace(x=greenEdges1$x, y=greenEdges1$y,z=greenEdges1$z, type="scatter3d", 
            mode="lines", name="GreenReflect", showscale=F,
            line = list(color=greenEdges1$line$color, width=greenEdges1$line$width)) %>%
  add_trace(x=redEdges$x, y=redEdges$y, z=redEdges$z, type="scatter3d", 
            mode="lines", name="RedLight",
            line=list(color=redEdges$line$color, width=redEdges$line$width), showscale=F) %>%
  add_trace(x=redEdges1$x, y=redEdges1$y, z=redEdges1$z, type="scatter3d", 
            mode="lines", name="RedReflect", showscale=F,
            line = list(color=redEdges1$line$color, width=redEdges1$line$width)) %>%
  # add a Ball centered at arg_min
  add_trace(x=1, y=1, z=0, type="scatter3d", mode="markers", 
            marker=list(size=10, color="blue"), name="point") %>%
  # add 2 extra balls and 2 extra tick reflection lines for G and B colors
  add_trace(x=greenEdges1$x[1], y=greenEdges1$y[1], z=0, type="scatter3d", 
            mode="markers", marker=list(size=10, color="green"), name="point", visible=F) %>%
  add_trace(x=redEdges1$x[1], y=redEdges1$y[1], z=0, type="scatter3d", 
            mode="markers", marker=list(size=10, color="red"), name="point", visible=F) %>%
  add_trace(x=greenEdges1$x, y=greenEdges1$y, z=greenEdges1$z, type="scatter3d", 
            mode="lines", name="GreenLight", showscale=F, visible=F,
            line = list(color=greenEdges1$line$color, width=blueEdges1$line$width)) %>%
  add_trace(x=redEdges1$x, y=redEdges1$y, z=redEdges1$z, type="scatter3d", 
            mode="lines", name="RedLight", showscale=F, visible=F,
            line = list(color=redEdges1$line$color, width=blueEdges1$line$width)) %>%
  # lastly add a dim light blue reflection line
  add_trace(x=blueEdges1$x, y=blueEdges1$y, z=blueEdges1$z, type="scatter3d", 
            mode="lines", name="BlueLight", showscale=F, visible=F,
            line = list(color=blueEdges1$line$color, width=greenEdges1$line$width)) %>%
  layout(legend = list(orientation = 'h'), 
         title="Object Surface Absorbs all Wavelengths except \n the ones Reflected out as Visual Colors",
         scene = list(aspectratio = list(x = 1,y = 1,z = 1)),
         updatemenus = updateColorMenus)

Let’s look at the electromagnetic wavelength spectrum. Solar radiation is radiant electromagnetic energy from the sun, which provides light, heat, and energy for photosynthesis. There are three important bands of solar radiation - ultraviolet, visible, and infrared. Infrared radiation makes up \(49.4\%\) of the energy reaching Earth surface, visible light accounts for another \(42.3\%\), and ultraviolet radiation makes up about \(8\%\) of the total solar radiation.

The amount and intensity of solar radiation impacting a specific location depends on latitude, altitude, season, time of day, cloud cover, and atmospheric pressure. Much of the radiation emitted from the sun is absorbed, reflected, or scattered in the atmosphere. Below we use data from National Renewable Energy Laboratory to show the Earth solar radiation.

wavelenthData <- read.csv("https://umich.instructure.com/files/26150288/download?download_frd=1", 
                          header = TRUE)    

fitNorm <- function(N, mean = 0, sd=1) {
  x <- rnorm(N, mean = mean, sd=sd) 
  density(x, bw = 30)   
}

normBlue <- fitNorm(10000, mean = 420, sd=70)
normGreen <- fitNorm(10000, mean = 550, sd=60)
normYellow <- fitNorm(10000, mean = 610, sd=70)

vline <- function(x = 0, color = "green") {
  list(type="line", y0=0, y1=1, yref="paper", x0=x, x1=x, line=list(color=color, dash="dot"))
}

plot_ly(data=wavelenthData, x=~Wavelength_nm, y=~OutsideAtm_EtrW.m.2.nm.1, 
        type="scatter", mode="lines", name="Outside Atmosphere") %>%
  add_trace(y=~Ground_DirectCircumsolar_W.m.2.nm.1, name="Sea Level") %>%
  # visible blue light density
  add_trace(x = normBlue$x, y = 50*normBlue$y, type = "scatter", mode = "lines", opacity=0.1,   
              line=list(color="black"), fill = "tozeroy", fillcolor="blue", name="VisibleBlue") %>%
  # visible green light density
  add_trace(x = normGreen$x, y = 80*normGreen$y, type = "scatter", mode = "lines", opacity=0.1, 
              line=list(color="black"), fill = "tozeroy", fillcolor="green", name="VisibleGreen") %>%
  # visible blue light density
  add_trace(x = normYellow$x, y = 60*normYellow$y, type = "scatter", mode = "lines", opacity=0.1,   
              line=list(color="black"), fill = "tozeroy", fillcolor="yellow", name="VisibleYellow") %>%
  layout(legend = list(orientation = 'h'), title="Solar Radiation Reaching Earth (visible and infrared light)",
         scene = list(aspectratio = list(x = 1,y = 1,z = 1)),
         xaxis=list(title="Wavelength (nm)"), 
         yaxis=list(title="Spectral iradiance intensity (W*(m^-2)*(nm^-1))"),
         shapes = list(vline(379), list(type="rect",text="Gamma, X, UV", #add visible spectrum
                                        fillcolor="red", # colorRamp(c("violet", "yellow", "red")),
                                        line = list(color = "red"), opacity = 0.2, y0 = 0.01, 
                                        y1 = 2.3, x0 = 380, x1 = 750),
                       vline(751), list(type="rect", text="Visible", fillcolor="gray", # add visible spectrum
                                        line = list(color = "gray"), opacity = 0.2, y0 = 0.01, 
                                        y1 = 2.3, x0 = 5, x1 = 370),
                       list(type = "rect", text="IR RW", fillcolor = "orange", # add visible spectrum
                                        line = list(color = "gray"), opacity = 0.2, y0 = 0.01, 
                                        y1 = 2.3, x0 = 390, x1 = 4000))) %>%
  add_text(showlegend = FALSE, x = 570, y = 0.6, text = "Visible Light") %>%
  add_annotations(showlegend = FALSE, x = 150, y = 0.1, text = "Gamma, X, UV",
           textfont = list(angle = 90, textangle = 90, orientation = 90, rotate = 90)) %>%
  add_annotations(showlegend = FALSE, x = 2500, y = 0.1, text = "IR, Radio waves",
           textfont = list(angle = 90, textangle = 90, orientation = 90, rotate = 90))

A simple theoretical model of color vision poses that the three primary RGB colors correspond to the three types of cones and all hues (RGB mixtures) can be untangled in the brain as linear combinations of stimulations arriving from the three types of cones. Computer screens and television rely on such three-color visual models. However, there is no unique set of three primary colors, e.g., yellow, green, and blue can also be used to decompose the visible spectrum into three independent components. This suggests that more complex color vision theoretical model may be useful to represent different combinations generating the same hue.

the figure above shows sensitivity curve models as y-axis intensities for given wavelength (\(\lambda\)) on the x-axis. The three model distributions corresponding to the retinal cones can be symmetric or skewed and cover specific wavelength ranges within the visible spectrum, \(350-700\ nm\). These models may be tuned by using evidence (observational data) for each of the three types of cones derived from direct measurements in animal and human eyes.

3 Forces and Equilibria

3.1 Forces

In 1D, a time-independent force \(F=F(x(t))\) moving a single particle of mass \(m\) in \(\mathbb{R}\) is conservative when it represents the gradient of some smooth real-valued function \(−V(x):\mathbb{R} \to \mathbb{R}\), i.e., \(F=\frac{d(-V)}{dx}\). This function \(V(x(t))\) is called the potential energy at time \(t\) and is determined only up to a real constant \(a\), as \(\frac{d(−V + a)}{dx}=-\frac{dV}{dx}\).

By Newton’s second law of motion, the particle traverses a path expressed as \[x(t):\mathbb{R}\to \mathbb{R},\] where the position function \(x(t)\in \mathbb{R}\) satisfies the following second order differential equation

\[m \frac{d^2x(t)}{dx^2} = −\frac{dV(x(t))}{dx}.\] The kinetic energy at time \(t\) is \[K(x(t))=\frac{1}{2}mv^2(t)\equiv \frac{1}{2}m \frac{dx(t)}{dt} \cdot \frac{dx(t)}{dt}\] Finally, the sum of the potential \(V\) and kinetic \(K\) energies defined the total energy at time \(t\), \[\underbrace{E(x(t))}_{total} = V(x(t)) + K(x(t))\equiv \underbrace{V(x(t))}_{potential}+ \underbrace{\frac{1}{2}m \left |\frac{dx(t)}{dt}\right |^2}_{kinetic}.\]

The force is conservative, since the total energy \(E(x(t))\) remains static (conservation of energy principle) despite the fact that as the time varies, both the kinetic \(K(x(t))\) and potential \(V(x(t))\) energies change along the path \(x(t)\):

\[\frac{dE(x(t))}{dt} = \frac{d}{dt}\left (\frac{1}{2}m \frac{dx(t)}{dt}\cdot \frac{dx(t)}{dt} +V(x(t)) \right)=\] \[2\left (\frac{1}{2}m \frac{dx(t)}{dt}\cdot \frac{d^2x(t)}{dt^2} \right)+\underbrace{\frac{d}{dt}(V(x(t)))}_{\frac{d}{dx}V(x) \cdot \frac{d}{dt}(x(t))}=\] \[\frac{d}{dt}(x(t)) \left ( \underbrace{m \frac{d^2x(t)}{dt^2} +\frac{d}{dx}V(x)}_{0} \right)=0.\]

When conservative forces are applied, the total energy \(E\) is constant (conserved) during the motion. The property that some physical systems may evolve in time but conserve some observable quantities is important for characterizing the systems.

Consider the example of having a mass attached to a spring \(x(t) \equiv q(t)\) and since \(F(x)=-kx=-\frac{d}{dx}\left (\frac{1}{2}kx^2 \right )\), the potential energy is \(V(q) = \frac{1}{2}kx^2\), and the kinetic energy is \(V(x)\). Hence, the total energy of the mass at time \(t\) is

\[E=\underbrace{\frac{1}{2}m\left |\frac{dx(t)}{dt}\right |^2}_{K} + \underbrace{\frac{1}{2}kx^2(t)}_{V}.\]

Under initial conditions \(V(0) = \frac{dV}{dx}(0) = 0\) and \(\frac{d^2V}{dx^2}(0) \gt 0\), the potential energy \(V\) attains a minimum, \(\min_{x}{V(x)}=V(0)=0\), which corresponds the force \(F(x) = −\frac{dV}{dx}\) that is trivial at \(x = 0\), i.e., the spring system is at a stable equilibrium point, \(x = 0\).

Let’s explore further this simple 1D spring system, subject to a (restoring) force \(F(x) = −\frac{dV}{dx}\), where the potential energy \(V(x)\) has a minimum at \(x = 0\), i.e., \(V(0) = \frac{dV}{dx}(0) = 0\) and \(\frac{d^2V}{dx^2}(0) \gt 0\).

Then, the Taylor series expansion for \(V(x)\) near at \(x = 0\) is \[V(x) = \underbrace{V(0)}_{0} + \underbrace{V'(0)}_{0}x + \frac{V''(0)}{2!}x^2+\cdots +\frac{V^{(n)}(0)}{n!}x^n+\cdots \ .\]

Within a small neighborhood around \(x=0\), we can approximate the potential energy by \(V(x)\approx \frac{1}{2}\ \underbrace{V''(0)}_{k\gt 0}\ x^2\), which is similar to the harmonic oscillator. So, locally around \(x=0\), the spring motion resembles a harmonic oscillator. Of course, the equilibrium point can be at any other location, e.g., \(x = x_o\), and the same problem formulation holds via Taylor series expansion near \(x = x_o\).

The potential function is determined only up to an additive constant, \(a\), which allows control over where \(V\) vanishes, e.g., \(V(x_o)=0\) and \(k=V''(x_o)\gt 0\) yields a stable system equilibrium near \(x_o\).

The importance of the general harmonic oscillator is realized as a model for any 1D system with a stable equilibrium point must behaves like a harmonic oscillator near its equilibrium point.

3.2 Translation

Rigid body translation (defined as a type of rigid transformation, since the Euclidean distance between every pair of points is preserved when applied) can be defined mathematically as a matrix addition between the origin and a translation vector, where \(\Delta x\), \(\Delta y\), and \(\Delta z\) are the displacements in each Cartesian direction:

\[\mathbf{t} = \left[\begin{array} {rrr} \Delta x & \Delta y & \Delta z \\ \end{array}\right].\]

3.3 Rotation

Rigid body rotations (also defined as a type of rigid transformation) can be defined mathematically as a matrix multiplication between the original matrix of points and a rotation matrix. Below are several examples of fundamental rotation matrices having a single parameter \(\theta\) specifying the angle of rotation around the given principal axis:

\[\mathbf{R_x} = \left[\begin{array} {rrr} 1 & 0 & 0 \\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \end{array}\right] , \]

\[\mathbf{R_y} = \left[\begin{array} {rrr} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{array}\right] , \]

\[\mathbf{R_z} = \left[\begin{array} {rrr} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{array}\right].\]

The following is a visualization of a translated and rotated spring using the linear algebra transformations as described above:

library(plotly)

#define translation vector in 3D with respect to the origin
delX <- 1
delY <- 2
delZ <- 3

t <- seq(-10, 10, 0.01)
curve <- data.frame(
  x = cos(10*t) + delX,
  y = sin(10*t) + delY,
  z = t + delZ
)

curve_mat <- as.matrix(curve)

# define rotation matrix in 3D and apply to data
theta_x <- pi/2 # radians, change this value to define other rotations
theta_y <- pi/4
theta_z <- pi/3

rot_matrix_x <- matrix(c(1, 0, 0, 0, cos(theta_x), sin(theta_x), 0, -sin(theta_x), cos(theta_x)), nrow=3)

rot_matrix_y <- matrix(c(cos(theta_y), 0, -sin(theta_y), 0, 1, 0, sin(theta_y), 0, cos(theta_y)), nrow=3)

rot_matrix_z <- matrix(c(cos(theta_z), sin(theta_z), 0, -sin(theta_z), cos(theta_z), 0, 0,0,1), nrow=3)

# apply transformation (rotation) to spring object
rotated_curveX <- curve_mat %*% rot_matrix_x
rotated_curveY <- curve_mat %*% rot_matrix_y 
rotated_curveZ <- curve_mat %*% rot_matrix_z 

rotated_curveZ_Y_X <- curve_mat %*% rot_matrix_x %*% rot_matrix_y %*% rot_matrix_z

std <- list(range = c(-4,4))

plot_ly(x = curve_mat[,1], y = curve_mat[,2], z = curve_mat[,3], type = "scatter3d", 
        mode = "lines", line=list(width=10), name="Original Spring") %>%
  add_trace(x = rotated_curveX[,1], y = rotated_curveX[,2], z = rotated_curveX[,3], 
            type = "scatter3d", mode = "lines", line=list(width=3), name="R_x") %>%
  add_trace(x = rotated_curveY[,1], y = rotated_curveY[,2], z = rotated_curveY[,3], 
            type = "scatter3d", mode = "lines", line=list(width=3), name="R_y") %>%
  add_trace(x = rotated_curveZ[,1], y = rotated_curveZ[,2], z = rotated_curveZ[,3], 
            type = "scatter3d", mode = "lines", line=list(width=3), name="R_z") %>%
  add_trace(x=rotated_curveZ_Y_X[,1], y=rotated_curveZ_Y_X[,2], z=rotated_curveZ_Y_X[,3], 
            type = "scatter3d", mode = "lines", line=list(width=3), name="R_z(R_y(R_x))") %>%
  layout(scene = list(xaxis = std, yaxis = std, zaxis = std, aspectmode = 'manual',
                      aspectratio = list(x = 1, y = 1, z = 1)), title="Spring Rotations",
         legend = list(orientation = 'h'))

4 Displacement, Velocity, and Acceleration

These concepts were already introduced in Chapter 1 (Math Foundations), including an example of a projectile motion model. Let’s now present one explicit example of sound waves.

4.1 Sound Waves

All of the five human senses track some form of waves. For instance hearing interprets sound waves as audible signals. Sound propagates through a medium (not through vacuum) as pressure waves displacing molecules in the medium, which can be gas, liquid, or solid. The simplest form of sound waves oscillate between compression and expansion of particle displacement in direction of wave propagation, parallel to the displacement amplitude. Suppose the horizontal x-axis represents the direction of wave propagation, the displacement amplitude, and the velocity molecules in the medium.

library(plotly)
size=300
t <- seq(from=0, to=2*pi, length.out=size)
pressure <- cos(5*t)
displacement <- sin(5*t)
y <- seq(from=-1, to=1, length.out=size)
z <- pressure %o% rep(1, size)
plot_ly(x=~t, y=~y, z = ~t(z), type = "heatmap", colors = "Greys", showscale= F) %>%
  add_trace(x=~t, y=~pressure, type="scatter", mode="lines", name="Pressure",
            line = list(color = 'blue', width = 5)) %>%
  add_trace(x=~t, y=~displacement, type="scatter", mode="lines", name="Displacement",
            line = list(color = 'red', width = 5)) %>%
  layout(title = "Displacement and Pressure of a Popagating Sound Wave", 
           xaxis=list(title="time"), yaxis=list(title="Value", scaleanchor="x", scaleratio=2.2),
           legend = list(title=list(text='<b>Functions</b>'), orientation = 'h')) # 1:1:1 aspect ratio
## Warning: 'scatter' objects don't have these attributes: 'z', 'showscale'
## Valid attributes include:
## 'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'fillpattern', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

## Warning: 'scatter' objects don't have these attributes: 'z', 'showscale'
## Valid attributes include:
## 'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'fillpattern', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

A harmonic wave function of the molecule position \(x\) and time \(t\) can be used to a model for the particle displacement \(\delta(x, t)= \delta_o \sin(\kappa x − \omega t)\), where the model parameters \(\kappa = \frac{2\pi}{\lambda}\) is the wave-number and \(\omega\) is the frequency of the sound wave propagating in the horizontal direction.

This harmonic wave function represents a solution of the 1D wave equation

\[{\frac {\partial^{2} \delta}{\partial t^{2}}} = c^{2} {\frac {\partial^{2} \delta}{\partial x^{2}}},\]

where \(c=\frac{\omega}{\kappa}\equiv \sqrt{\frac{E}{\rho}}\) is the phase velocity (speed) of the sound wave, \(E\) is Young’s modulus (compression modulus), \(rho\) is the density of the ambient material. Let’s consider two examples sound propagation in gas and liquid.

This unidirectional wave equation model of sound propagation generalizes to 4D Minkowski spacetime as vector wave equation, which naturally reduces to its scalar wave equation counterpart in a single spatial dimension. Assume the medium is continuous and homogeneous, \(\boldsymbol {x}\in \mathbb{R}^3\) is the space coordinate vector, and \(E\) is the constant Young’s modulus of compression elasticity. Then, the elastic deflection \(\delta( \boldsymbol {x} , t )\) is associated with the stress tensor \(\boldsymbol {T}=E\nabla {\delta}\) whose local equilibrium is the tension force \(div{\boldsymbol {T} = \nabla \cdot (E\nabla \delta)=E\Delta {\boldsymbol {\delta}}}\), where \(\nabla\) and \(\Delta\) are the the (first-order) gradient and the (second-order) Laplacian operators.

Therefore, the wave equation can be expressed as equilibrium equating the inertial force \(\rho \frac{\partial^{2}{\boldsymbol {u}}}{\partial t^{2}}\) due to the local acceleration in time \(\frac{\partial^{2}{\boldsymbol {u}}}{\partial t^{2}}\) and the tension force \(div{\boldsymbol {T} = \nabla \cdot (E\nabla \delta)=E\Delta {\boldsymbol {\delta}}}\) (divergence) associated with the stress tensor \(T\).

\[\underbrace{\rho \frac{\partial^{2}\delta( \boldsymbol {x} , t )} {\partial t^{2}}}_{inertial\ force} = \underbrace{E\Delta_{\boldsymbol {x}} {\delta( \boldsymbol {x} , t )}}_{tension\ force}.\] Simplifying the equation using the phase velocity of the sound wave, \(c=\frac{\omega}{\kappa}\equiv \sqrt{\frac{E}{\rho}}\), we get

\[{\frac {\partial ^{2}{\delta}}{\partial t^{2}}}-c^{2}\Delta {\delta}=\underbrace{\boldsymbol {0}}_{vector}.\]

This second-order vector PDE has a pair of independent solutions corresponding to two waves travelling in opposite directions \(\pm c\), since the quadratic phase velocity term \(c^{2}=(+c)^{2}=(-c)^{2}\) yields two opposite momenta.

Recall that adiabatic and isothermal processes represent to types of thermodynamic processes. In adiabatic processes there is no heat supply to or release from the body which undergoes a thermodynamic state change, i.e., the system is in adiabatic isolation; however, the system temperature can vary to accommodate the lack of any heat transfer. In isothermal processes, the temperature of the body remains fixed during the thermodynamic process, however, other system parameters (e.g., volume, pressure) may have process-specific changes.

In an ideal gas medium, the compression modulus can be expressed as \(E=p\) (for isothermic processes) or \(E = \gamma p\) (for adiabatic processes), where \(\gamma=\frac{\text{heat at constant pressure}}{\text{heat at constant volume}}\) and \(p\) is the ideal gas pressure.

For example, sound waves propagating in air are modeled as adiabatic processes since the thermal equilibrium by heat exchange cannot be reached for any frequency above \(20 Hz\) (normal human hearing range is \(20Hz-20kHz\)). To summarize, the sound wave velocity in gas and liquid media are

\[c_{gas}=\sqrt{\frac{\gamma p}{\rho}} \ \ \ \& \ \ \ c_{liquid}=\sqrt{\frac{E}{\rho}}.\]

At normal air conditions of pressure and humidity, the sound speed is \(c_{air} \sim 330 m/s\), whereas in water (a much denser medium), the sound speed is about five times faster, \(c_{water} \sim 1500 m/s\).

Hooke’s law suggests that the force \(F\) needed to extend or compress a spring by some distance \(x\) scales linearly with respect to that distance, i.e., \(F_s = \kappa x\), where the linear scaling parameter \(\kappa\) is a constant dependent on the spring characteristics, e.g., its stiffness, and \(x\) is small compared to the total possible deformation of the spring.

Let’s explicate the sound wave amplitude velocity of molecules using the first-order derivative of the displacement, \(\delta(x, t)\), \(u_x = \frac{\partial \delta (x, t)}{\partial t} = \underbrace{−\omega \delta_o}_{u_o} \cos(\kappa x − \omega t)\). Locally, Hooke’s law describes the propagation of the pressure wave \(p_x = −E \frac{\partial \delta (x, t)}{\partial x} = \underbrace{−E \delta_o \kappa}_{p_o} \cos(kx − \omega t)\), where \(u_o\) and \(p_o\) are the velocity and pressure amplitudes, respectively.

For a given area of sound propagation, \(A\), the expected power \(P\) of the sound wave and its expected intensity \(I\) are computed as time averages of the inner-product of the force \(F_x\) and the particle velocity \(u_x\)

\[\underbrace{\langle P\rangle}_{mean \ power} = \langle F_x u_x \rangle \ \ \ \& \ \ \ \underbrace{\langle I\rangle}_{mean \ intensity}\equiv \left\langle \frac{P}{A}\right\rangle = \left\langle \frac{F_x}{A} u_x \right\rangle = \langle p_x u_x \rangle = E\omega \kappa \delta_o \underbrace{\langle \cos^2(\kappa x - \omega t) \rangle}_{\frac{1}{2}}= \frac{1}{2} E\omega \kappa \delta_o .\]

The last equality reflects the fact that over one period \(T=\frac{2\pi}{\omega}\) for a fixed position \(x\),
\(\int_{0}^{T}\cos^2(\kappa x - \omega t)dt \equiv \int_{0}^{T}\sin^2(\kappa x - \omega t)dt=\frac{1}{2}\). This is derived as follows. Since \(\cos(2x) = \cos^2(x) - \sin^2(x)\) and \(\sin^2(x) + \cos^2(x) = 1\), we conclude that \(\cos(2x) = 2\cos^2(x) -1\) and \(\cos^2(x)=\frac{\cos(2x)+1}{2}\). We can use the reverse chain rule to compute

\[\int_{0}^{T}\cos^2(\kappa x - \omega t)dt = \frac{1}{2}\int_{0}^{T}(\cos(2(\kappa x - \omega t))+1)dt = \frac{1}{2} \left (t - \frac{\sin(2(\kappa x - \omega t))}{\omega}\right ) \Big |_{0}^{T} = \frac{T}{2} .\]

CHECK: Why are we off by a factor of the period \(T\)???

The sound intensity \(\langle I \rangle=\frac{1}{2} E\omega \kappa \delta_o \equiv \frac{1}{2} (\rho c) u_o^2\) is measured in \(Watts\times m^{-2}\) units. The acoustic impedance of the sound wave, \(\rho c\) is important for evaluating the sound transmission and reflection coefficients as the sound wave traverses between different media. Similarly to the refractive index in optics, the acoustic impedance characterizes the sound intensity in different media. Approximations of sound velocity, tissue density, and acoustic impedance for different body tissues are shown below.

Tissue Sound velocity (m s-1) Density (kg m-3) Acoustic impedance (kg m-2 s-1)
Air 330 1.3 430
Lung 350 240 0.18×106
Water 1500 998 1.5×106
Blood 1570 1060 1.65×106
Fat 1470 970 1.38×106
Brain 1540 1030 1.58×106
Muscle 1580 1070 1.71×106
Skin 1600 1100 1.62×106
Bone 3600 1700 7.8×106

The following 3D scene illustrates the direct proportionality of the sound wave intensity \(I\) to the change in the pressure squared and its inverse proportionality to the density and the speed. Suppose we start with an initial (undisturbed by sound) volume, which is then subjected to a sound wave at time \(t\).

The initial volume of the undisturbed 3D cube (blue ambient medium) is \(V=A\Delta x\), where \(A\) is the 2D base area of the cube, and \(\Delta x\) is the horizontal extent of the cube/parallelepiped. As the sound wave moves with time through the cube, volume is displaced and expanded/contracted (orange block). The volume change is \(\Delta V = A \Delta s =A(s_2−s_1)\), where \(s_1\) and \(s_2\) are the displacements of the block’s leading and trailing edges, respectively. In the 3D scene, \(s_2 > s_1\) indicating that cube is expanded by the wave at the given time, \(t\). Of course, the orange block may also be compressed by the sound wave, i.e., \(s_2 < s_1\), or be unchanged, i.e., \(s_2 \equiv s_1\).

The raw and fractional volume changes are \[\Delta V = \overbrace{A}^{\ \ \ Cube\\ face\ area} \Delta s =A(s_2−s_1)=A\left (s(\overbrace{x+\Delta x}^{Cube\ depth},t)−s(x,t)\right ),\] \[\frac{d V}{V} = \lim_{\Delta x\to 0} {\frac{A (s(x+\Delta x,t)−s(x,t))}{A\Delta x}}= \frac{\partial s(x,t)}{\partial x} .\]

The fractional volume change \(\frac{d V}{V}\) is related to the pressure and bulk modulus. The volume stress \(\Delta p =\frac{F}{A}\) is the ratio of the magnitude of the normal force \(F\) to the area \(A\). Whereas, the volume strain is the fractional change in volume, \(\frac{\Delta V}{V}\), and the bulk modulus \(B=-\frac{\Delta p}{\frac{\Delta V}{V}}\) is the ratio of the volume stress to the volume strain:

The negative sign guarantees that \(B >0\), i.e., forces pressing from the outside always cause the volume to decrease (compression), so that the fractional change in volume is negative. In other words, the minus sign reflects the inverse-proportionality of volume and pressure.

Sound propagation with time is often modeled as a pressure wave. Let’s look at the change in pressure from average pressure, \(\Delta p= \Delta p_{max} \sin(\kappa x\mp \omega t)\).

The maximum pressure change is denoted by \(\Delta p_{max}\), \(\kappa =\frac{2\pi}{\lambda}\) is the wave number, \(\omega=\frac{2\pi}{T}=2\pi f\) is the angular frequency, \(T\) is the period, and we assume a trivial initial phase \(\phi=0\). The wave speed is \(v=\frac{\omega}{\kappa}\equiv \frac{\lambda}{T}\).

Alternatively, sound waves can be modeled as displacements of molecules in the ambient environment. For instance, displacement of air molecules can be modeled using a cosine function \(s(x,t)=s_{max}\cos(\kappa x\mp \omega t)\). Then the pressure change caused by the sound wave is \(\Delta p(x,t) =−B \frac{\Delta V}{V}=−B\frac{\partial s(x,t)}{\partial x}=\) \(B\kappa s_{max}\sin(\kappa x - \omega t)=\Delta p_{max}\sin(\kappa x - \omega t)\).

The sound wave power is the product of force and the velocity of the wave oscillations of the medium as temporal rate of change of the displacement \(v(x,t)=\frac{\partial}{\partial t}s(x,t)=\frac{\partial}{\partial t}(s_{max}\cos(\kappa x−\omega t))=\omega s_{max}\sin(\kappa x−\omega t)\).

Hence, the wave intensity is \[I=\Delta p(x,t)v(x,t)= B\kappa s_{max}\sin((\kappa x−\omega t)[s_{max}\omega \sin((\kappa x−\omega t)]= B\kappa \omega s^2_{max}\sin^2(\kappa x−\omega t).\]

At a given spatial position \(x\) in the medium, the wave intensity averaged over one period \(T=\frac{2\pi}{\omega}\) is obtained by integrating the intensity over a single period, i.e., \(I=\frac{B\kappa \omega s_{max}^2}{2}\).

Finally, substituting \(\Delta p_{max}=B \kappa s_{max}\), \(c=\sqrt{\frac{B}{\rho}}\), and \(c=\frac{\omega}{\kappa}\), \[I=\frac{B\kappa \omega s_{max}^2}{2}=\frac{B^2\kappa^2 \omega s_{max}^2}{2B\kappa}= \frac{\omega (\Delta p_{max})^2}{2(\rho c^2)\kappa}=\frac{c (\Delta p_{max})^2}{2(\rho c^2)}= \frac{(\Delta p_{max})^2}{2\rho c}.\]

### Draw the 3D scene
points <- data.frame(cbind(x=c(0, 0, 1, 1, 0, 0, 1, 1), y=c(0, 1, 1, 0, 0, 1, 1, 0), z=c(0, 0, 0, 0, 1, 1, 1, 1)))
faces <- data.frame(cbind(i = c(7, 0, 0, 0, 4, 4, 6, 6, 4, 0, 3, 2),
                          j = c(3, 4, 1, 2, 5, 6, 5, 1, 0, 1, 7, 3),
                          k = c(0, 7, 2, 3, 6, 7, 1, 2, 5, 5, 6, 6)))
## i,j,k face indices
## vectors of vertex indices, i.e. integer values between 0 and the length of the vertex vectors, representing the 2-cell vertex of the triangle face. For example, {i[m], j[m], k[m]} together represent face m (triangle m) in the mesh, where j[m] = n points to the triplet {x[n], y[n], z[n]} in the vertex arrays. Therefore, each element in j represents a point in space, which is the second vertex of a triangle.

## Mock up Sound-wave front as a sphere (x-sphereXcenter)^2 + (y-0.5)^2 + (z-0.5)^2 = 2^2
# Define Half-Sphere Surface propagating in x-direction
R <- 0.5  #radius
rSeq = rep(R,100) # r = radius
u <- seq(from=0, to=pi/2, length.out = 100)
v <- seq(from=0, to=2*pi, length.out = 100)
z = R * cos(v) %o% sin(u) +0.5     # z = r2 cos(phi)sin(psi)
y = R * sin(v) %o% sin(u) +0.5     # y = r2 sin(phi)sin(psi)
x = R * (rSeq %o% cos(u)) -0.5
# plot_ly(x = ~x, y = ~y, z = ~z, type = 'surface', colors="gray", name="Sound-Front",
#         hoverinfo="Sound-Front", legendshow=F, opacity=0.6) %>%
#   add_trace(x = x+0.5, y = ~y, z = ~z, type = 'surface', colors="green", name="Wave-Propagator",
#         hoverinfo="Sound-Front", legendshow=F, opacity=0.6)

  
plot_ly(type = 'mesh3d', opacity=0.3, name="Original Volume", hovertemplate = paste('(%{x}, %{y}, %{z})'),
  # Undisturbed 3D Cube
        x = points[,1], y = points[,2], z = points[,3],  # vertices
        i = faces[,1], j = faces[,2], k = faces[,3]) %>% # i, j and k give the vertices of triangles
  add_trace(type = 'mesh3d', opacity=0.5, name="2D Base Plane Area", hovertemplate = paste('(%{x}, %{y}, %{z})'),
  # 2D Base Area: faces 10=(0,1,5) and 9=(0,4,5)
        x = points[,1], y = points[,2], z = points[,3], 
        i = faces[c(9,10),1], j = faces[c(9,10),2], k = faces[c(9,10),3]) %>%
  add_trace(type = 'mesh3d', opacity=0.5, name="Sound-displaced/expanded Volume", 
            hovertemplate = paste('(%{x}, %{y}, %{z})'), 
        # vertices
        x = (1.3)*(points[,1]+0.6), y = points[,2], z = points[,3], # expansion-factor=1.3, displacement=0.6
        # i, j and k give the vertices of triangles
        i = faces[,1], j = faces[,2], k = faces[,3]) %>% 
  # Add the sound-wave front propagator: 1
  add_trace(x = x+0.2, y = ~y, z = ~z, type = 'surface', colors="gray", name="Wave-Propagator 1",
            hoverinfo="Sound-Front 1", showlegend=F, opacity=0.1, showscale = FALSE) %>% 
  # Add the sound-wave front propagator
  add_trace(x = x+0.4, y = ~y, z = ~z, type = 'surface', colors="gray", name="Wave-Propagator",
            hoverinfo="Sound-Front", showlegend=F, opacity=0.6, showscale = FALSE) %>%
  # Add the sound-wave front propagator
  add_trace(x = x+0.6, y = ~y, z = ~z, type = 'surface', colors="gray", name="Wave-Propagator 2",
            hoverinfo="Sound-Front 2", showlegend=F, opacity=0.1, showscale = FALSE) %>%
  # Add s1, the displacement of the Leading edge of the Cube
  add_trace(x=c(0,(1.3)*(points[1,1]+0.6)), y=0, z=0, type="scatter3d", mode="lines+markers",opacity=1, 
            line = list(width = 20, color="blue"), name="s1 Leading-Plane Displacement", showlegend=F) %>% 
  # Add s2, the displacement of the Trailing edge of the Cube
  add_trace(x=c(1,(1.3)*(points[length(points$x),1]+0.6)), y=0, z=0, type="scatter3d", mode="lines+markers", opacity=1, 
            line = list(width = 20, color="orange"), name="s2 Trailing-Plane Displacement", showlegend=F) %>% 
  # Add Delta s=s2-s1 difference
  add_trace(x=c((1.3)*(points[1,1]+0.6),1), y=0.5, z=1, type="scatter3d", mode="lines+markers",opacity=1, 
            line = list(width = 20, color="red"), name="Ds=s2-s1 Difference", showlegend=F) %>% 
  # Add Delta x
  add_trace(x = c(0,1), y = c(0,0), z = c(0.2,0.2), type = "scatter3d", mode = "lines_markers", 
            line = list(width = 10, color="yellow"), name = "Delta x", showlegend = FALSE, opacity=1) %>%
  layout(scene = list(aspectratio = list(x = 2, y = 1, z = 1), dragmode = "turntable", showscale = FALSE,
      annotations = list(list(showarrow = F, x = 0.3, y = 0.5, z = 0.5, text = "Undisturbed Volume",
                              xanchor = "middle", xshift = 10, opacity = 1.0), 
                         list(showarrow = F, x = 1.6, y = 0.5, z = 0.5, text = "Sound-disturbed Volume",
                              xanchor = "middle", xshift = 10, opacity = 1.0), 
                         list( x = 0.0, y = 0.5, z = 0.5, text = "(Base Plane) Area",
                              textangle = 0, ax = 0, ay = -75, 
                              font = list(color = "black", size = 12),
                              arrowcolor = "black", arrowsize = 3, arrowwidth = 1,arrowhead = 1),
                         list( x = 0.5, y = 0.0, z = 0.2, text = "Delta x (Dx)",
                              textangle = 0, ax = 0, ay = -75, 
                              font = list(color = "black", size = 12),
                              arrowcolor = "black", arrowsize = 3, arrowwidth = 1,arrowhead = 1),
                         list( x = 0.85, y = 0.5, z = 1.0, text = "Displacement Change Delta (Ds=s2-s1)",
                              textangle = 0, ax = 0, ay = -75, 
                              font = list(color = "black", size = 12),
                              arrowcolor = "black", arrowsize = 3, arrowwidth = 1,arrowhead = 1)),
      camera = list(eye = list(x = -1.5, y = -1.5, z = 1.5))))
## Warning: 'surface' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## '_deprecated', 'autocolorscale', 'cauto', 'cmax', 'cmid', 'cmin', 'coloraxis', 'colorbar', 'colorscale', 'connectgaps', 'contours', 'customdata', 'customdatasrc', 'hidesurface', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'lighting', 'lightposition', 'meta', 'metasrc', 'name', 'opacity', 'opacityscale', 'reversescale', 'scene', 'showlegend', 'showscale', 'stream', 'surfacecolor', 'surfacecolorsrc', 'text', 'textsrc', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

## Warning: 'surface' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## '_deprecated', 'autocolorscale', 'cauto', 'cmax', 'cmid', 'cmin', 'coloraxis', 'colorbar', 'colorscale', 'connectgaps', 'contours', 'customdata', 'customdatasrc', 'hidesurface', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'lighting', 'lightposition', 'meta', 'metasrc', 'name', 'opacity', 'opacityscale', 'reversescale', 'scene', 'showlegend', 'showscale', 'stream', 'surfacecolor', 'surfacecolorsrc', 'text', 'textsrc', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

## Warning: 'surface' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## '_deprecated', 'autocolorscale', 'cauto', 'cmax', 'cmid', 'cmin', 'coloraxis', 'colorbar', 'colorscale', 'connectgaps', 'contours', 'customdata', 'customdatasrc', 'hidesurface', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'lighting', 'lightposition', 'meta', 'metasrc', 'name', 'opacity', 'opacityscale', 'reversescale', 'scene', 'showlegend', 'showscale', 'stream', 'surfacecolor', 'surfacecolorsrc', 'text', 'textsrc', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'
## Warning: 'scatter3d' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## 'connectgaps', 'customdata', 'customdatasrc', 'error_x', 'error_y', 'error_z', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'projection', 'scene', 'showlegend', 'stream', 'surfaceaxis', 'surfacecolor', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

## Warning: 'scatter3d' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## 'connectgaps', 'customdata', 'customdatasrc', 'error_x', 'error_y', 'error_z', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'projection', 'scene', 'showlegend', 'stream', 'surfaceaxis', 'surfacecolor', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

## Warning: 'scatter3d' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## 'connectgaps', 'customdata', 'customdatasrc', 'error_x', 'error_y', 'error_z', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'projection', 'scene', 'showlegend', 'stream', 'surfaceaxis', 'surfacecolor', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

## Warning: 'scatter3d' objects don't have these attributes: 'i', 'j', 'k'
## Valid attributes include:
## 'connectgaps', 'customdata', 'customdatasrc', 'error_x', 'error_y', 'error_z', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'projection', 'scene', 'showlegend', 'stream', 'surfaceaxis', 'surfacecolor', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

4.2 Hearing

The temporal bone in mammalian skull contains the smallest elements of the muscular-skeleton which control balance, posture, stable gaze, and hearing. The complex evolution of the jaw joint into the contemporary mammalian ear ossicles is not fully known. The ear anatomy in birds and reptiles is much simpler than the mammalian hearing apparatus.

The outer ear canal channels sound to the eardrum and is responsible for the sound sensitivity in the \(2-5kHz\) range. The middle ear transforms sound-waves into mechanical vibrations on the cochlea, which further converts these signals into electrical impulses communicated via nerve synaptic connections to the auditory cortex in the brain temporal lobe.

The sensitivity of the human ear to sound intensities is sound-frequency dependent. In general, at frequency \(1 kHz\), the minimum (detectable sound) and maximum (pain or hearing loss) of the sound intensity are \(I_o=10^{−12}W/m^2\) and \(I_{max}=1W/m^2\). Using the log transformation of the normalized sound intensity value allows us to compress the wide range of intensities into something we can interpret more easily. Often, sound intensity is reported in decibels, \(\beta (dB)\), instead of Watts per meter squared (\(W/m^2\)), where deci refers to the decimal and bel is in honor of Alexander Graham Bell:

\[\beta (dB)=10\log_{10}\left (\frac{I}{I_o}\right ).\]

The decibel measure of sound intensity is a fraction, and therefore, is a unitless quantity, relativized to a fixed \(I_o=10^{−12}W/m^2\) sound level. Despite the fact that air molecules disrupted by a sound wave of intensity \(I_o\) vibrate only over a short distance (less than one molecular diameter, \(10^{-8}cm\)) and the gauge pressure of these vibrationsare minute (\(<10^{-9}atm\)), normal human hearing is still sensitive to minimal sound intensities around \(10^{-12} W/m^2\). The summary table below explicates the spectrum of sound intensities and their impact.

β (dB) Intensity I (W/m2) Example/effect
0 1×10-12 Hearing lower-bound at 1000 Hz
10 1×10-11 Normal breathing
20 1×10-10 Whisper at 1m distance
30 1×10-9 Whispering 0.1m distance
40 1×10-8 Home background
50 1×10-7 Office background
60 1×10-6 Directed Normal conversation
70 1×10-5 Noisy office
80 1×10-4 Loud home appliences
90 1×10-3 Train, airplane
100 1×10-2 Noisy factory, sirens at 50m, hearing loss after 10+ hours of exposure
110 1×10-1 Car horn, hearing loss from repetaed 30-min daily exposures
120 1 Loud rock concert; pneumatic chipper at 2m, initial pain
140 1×102 Jet airplane at 30m; severe pain, hearing loss in 1+ min
160 1×104 Damage to eardrums

4.2.1 Sound Intensity Example

Let’s try to estimate the sound intensity level (in decibels) for a sound wave with pressure amplitude of \(\Delta p=0.7 Pa\) propagating in normal air at \(0^o C\).

We’ll use the intensity-pressure equation \(I=\frac{(\Delta p)^2}{2\rho c}\) and the intensity-decibel relation \(\beta (dB)=10\log_{10}\left (\frac{I}{I_o}\right )\).

At normal sea-level atmospheric pressure and \(0^o C\) temperature, normal air density is \(1.29kg/m^3\) and sound waves travel at \(331 m/s\). Hence, \(I=\frac{(\Delta p)^2}{2\rho c}=\frac{(0.7Pa)^2}{2(1.29kg/m^3) 331 m/s}=5.737839\times 10^{-4} W/m^2\).

The corresponding decibel level is
\(\beta (dB)=10\log_{10}\left (\frac{I}{I_o}\right )=10\log_{10}\left (\frac{5.737839\times 10^{-4}}{10^{-12}}\right )=\) \(10\log_{10}(5.737839\times 10^{8} )=87.58748\). In logarithmic terms, this high level of sound intensity \(87.6dB\) corresponds to six times greater sound intensity of another \(80dB\) sound. A significant sound intensity amplification by a factor of six reflects only a marginal difference of \(7.6dB\) in the corresponding sound decibel level.

4.3 Allen-Neely Model

The Allen & Neely sound-intensity and frequency model can be used to visualize a sound intensity-frequency diagram showing the relations between audio frequencies and the sound intensity. the level-sets in this 2D plot expose the spectra of different human daily activities and relations between audibility, threshold of hearing pain, ranges for speech, music, and other audio perceptions.

This is a \(3^{rd}\) order regression model fit on partial data reported by Fletcher & Munson (Loudness, its definition, measurement, and calculation) in 1933.

The model coefficients reported in the table specify cubic polynomials for different frequency ranges.

\[\beta_{dB}(x|f) = c_o(f) + c_1(f) I + c_2(f) I^2 +c_3(f) I^3,\]

where \(f\) is the sound wave frequency in \(kHz\) and the model predicted value \(\beta\) is the sound level (in \(dB\)) at the specified sound intensity \(I\). At the lowest frequency range (\(62-125Hz\)), the model used cubic polynomials, whereas quadratic polynomials suffice for modeling the relation in the higher frequency range.

tabl <- "
| f(kHz) |  c<sub>3</sub>  |    c<sub>2</sub>  |c<sub>1</sub>  |    c<sub>0</sub>  |
|--------|:----------:|----------:|:----------:|----------:|
| 0.062 | 7.46169e-05 | -0.00984189 | 0.74629 | 0.425879 |
| 0.125 | 5.2594e-05 | -0.00654132 | 0.800557 | 0.295663 |
| 0.25 | 0 | 0.00124457 | 0.720323 | 0.780066 |
| 0.5 | 0 | 0.00209933 | 0.761911 | 0.467849 |
| 1 | 0 | 0 | 1 | 0 |
| 2 | 0 | -0.0011956 | 1.14141 | -0.622967 |
| 4 | 0 | -0.00240718 | 1.2314 | -0.393083 |
| 5.65 | 0 | -0.00272458 | 1.24014 | 0.481007 |
| 8 | 0 | -0.00232339 | 1.20659 | 0.0426691 |
| 11.3 | 0 | -0.002439 | 1.24474 | -1.51871 |
| 16 | 0 | -0.000566296 | 1.03446 | -1.82771 |
"
cat(tabl) # output the table in a format good for HTML/PDF/docx conversion
## 
## | f(kHz) |   c<sub>3</sub>  |    c<sub>2</sub>  |c<sub>1</sub>  |    c<sub>0</sub>  |
## |--------|:----------:|----------:|:----------:|----------:|
## | 0.062 | 7.46169e-05 | -0.00984189 | 0.74629 | 0.425879 |
## | 0.125 | 5.2594e-05 | -0.00654132 | 0.800557 | 0.295663 |
## | 0.25 | 0 | 0.00124457 | 0.720323 | 0.780066 |
## | 0.5 | 0 | 0.00209933 | 0.761911 | 0.467849 |
## | 1 | 0 | 0 | 1 | 0 |
## | 2 | 0 | -0.0011956 | 1.14141 | -0.622967 |
## | 4 | 0 | -0.00240718 | 1.2314 | -0.393083 |
## | 5.65 | 0 | -0.00272458 | 1.24014 | 0.481007 |
## | 8 | 0 | -0.00232339 | 1.20659 | 0.0426691 |
## | 11.3 | 0 | -0.002439 | 1.24474 | -1.51871 |
## | 16 | 0 | -0.000566296 | 1.03446 | -1.82771 |
modelCoef <- c(0.062, 7.46169e-05, -0.00984189, 0.74629, 0.425879, 
         0.125, 5.2594e-05, -0.00654132, 0.800557, 0.295663, 
         0.25, 0, 0.00124457, 0.720323, 0.780066, 
         0.5, 0, 0.00209933, 0.761911, 0.467849, 
         1, 0, 0, 1, 0, 
         2, 0, -0.0011956, 1.14141, -0.622967, 
         4, 0, -0.00240718, 1.2314, -0.393083, 
         5.65, 0, -0.00272458, 1.24014, 0.481007, 
         8, 0, -0.00232339, 1.20659, 0.0426691, 
         11.3, 0, -0.002439, 1.24474, -1.51871, 
         16, 0, -0.000566296, 1.03446, -1.82771)
mat <- t(matrix(modelCoef, nrow = 5, ncol = 11)) # transpose to get correct coeff table 
colnames(mat) <- c("f", "c3", "c2", "c1", "c0")
mat
##            f          c3           c2       c1         c0
##  [1,]  0.062 7.46169e-05 -0.009841890 0.746290  0.4258790
##  [2,]  0.125 5.25940e-05 -0.006541320 0.800557  0.2956630
##  [3,]  0.250 0.00000e+00  0.001244570 0.720323  0.7800660
##  [4,]  0.500 0.00000e+00  0.002099330 0.761911  0.4678490
##  [5,]  1.000 0.00000e+00  0.000000000 1.000000  0.0000000
##  [6,]  2.000 0.00000e+00 -0.001195600 1.141410 -0.6229670
##  [7,]  4.000 0.00000e+00 -0.002407180 1.231400 -0.3930830
##  [8,]  5.650 0.00000e+00 -0.002724580 1.240140  0.4810070
##  [9,]  8.000 0.00000e+00 -0.002323390 1.206590  0.0426691
## [10,] 11.300 0.00000e+00 -0.002439000 1.244740 -1.5187100
## [11,] 16.000 0.00000e+00 -0.000566296 1.034460 -1.8277100
# compute the frequency-intensity distribution (z) using the AN-model
freq_kHz <- seq(from=mat[1,1], to=mat[dim(mat)[1],1], length.out=200) # frequency uniform grid
soundIntensity_dB <- seq(from=0, to=130, length.out=200) # intensity uniform grid

z_fun <- function(f, I) { 
  z = 0
  # for each of the 11 AN-model frequency bands use appropriate poly-model
  for (i in 1:(dim(mat)[1]-1)) { 
    # this is the right model freq-range to get the poly coeffs (mat[,i])
    if (f[i] >= mat[i,1] & f[i] < mat[i+1,1]) {
      # \beta_{dB}(x|f) = c_o(f) + c_1(f) I + c_2(f) I^2 +c_3(f) I^3
      # z = mat[i,5] + mat[i,4]*I + mat[i,3]*(I^2) + mat[i,2]*(I^3) 
      z = mat[i,5] + mat[i,4]*I + mat[i,3]*(I^2) + mat[i,2]*(I^3)
    }
  }
  return (z)
}
# z_fun(freq_kHz[1], soundIntensity_dB[1]);z_fun(freq_kHz[1], soundIntensity_dB[1+1]);z_fun(freq_kHz[1+1], soundIntensity_dB[1]);
z <- outer(freq_kHz, soundIntensity_dB, FUN="z_fun")

z <- matrix(0, nrow = length(freq_kHz), ncol = length(soundIntensity_dB))

for (f in 1:dim(z)[1]) { 
  for (i in 1:(dim(mat)[1]-1))  {
    if (f>=mat[i,1] & f<=mat[i+1,1]) { # this is the right model freq-range to get the poly coeffs (mat[,i])
        # \beta_{dB}(x|f) = c_o(f) + c_1(f) I + c_2(f) I^2 +c_3(f) I^3
        # z = mat[i,5] + mat[i,4]*I + mat[i,3]*(I^2) + mat[i,2]*(I^3) 
      for (s in 1:dim(z)[2])  {
        z[f,s] = mat[i,5] + mat[i,4]*soundIntensity_dB[s] + 
                 mat[i,3]*(soundIntensity_dB[s]^2) + mat[i,2]*(soundIntensity_dB[s]^3)
      }
    }
  }
}

plot_ly() %>% add_trace(x=freq_kHz, y=soundIntensity_dB, z=t(z), type="surface", opacity=0.7) %>%
  # # add vertical axis at the minimum
  # add_trace(x=7, y=-8, z=c(-24:350), type="scatter3d", mode="lines", 
  #             line = list(width = 6, color="navy blue"), name="Z",
  #             hoverinfo="none") %>%
  # # add a Ball centered at arg_min
  # add_trace(x=7, y=-8, z=-24, type="scatter3d", mode="markers")
  layout(title="Allen & Neely sound-intensity and frequency model",
         scene = list(xaxis=list(title="Frequency (kHz), log-scale", tickmode = "array", nticks = dim(mat)[1],
                                 tickvals = mat[,1], type = "log"),
                      yaxis=list(title="Sound Level (dB)"), 
                      zaxis=list(title="AN Model")))
plot_ly() %>% add_trace(x=freq_kHz, y=soundIntensity_dB, z=t(z), type="heatmap", opacity=0.9) %>%
  layout(title="Allen & Neely sound-intensity and frequency model",
         xaxis=list(title="Frequency (kHz), log-scale", tickmode = "array", nticks = dim(mat)[1],
                    tickvals = mat[,1], type = "log"),
         yaxis=list(title="Sound Level (dB)"))
# this works!!!
m <- matrix(0, nrow = dim(mat)[1], ncol = length(soundIntensity_dB))
for (i in 1:(dim(mat)[1]))  {
   for (s in 1:length(soundIntensity_dB))  {
        m[i,s] = mat[i,5] + mat[i,4]*soundIntensity_dB[s] + 
                 mat[i,3]*(soundIntensity_dB[s]^2) + mat[i,2]*(soundIntensity_dB[s]^3)
  }
} # ;image(m)

# Define the parametric curve minimum auditory field (MAF)
xAudibleCurve <- c(mat[,1], rev(mat[1:length(mat[,1])-1,1]))
# correct left- and right-most points to keep curve in graph FOV
xAudibleCurve[1] <- xAudibleCurve[1]+0.01
xAudibleCurve[11] <- xAudibleCurve[11]-2; xAudibleCurve[21] <- xAudibleCurve[21]+0.01 
yAudibleCurve <- c(120, 70, 50, 30, 15, 5, -2, 0, 18, 40, 120, 120, 110, 110, 120, 120, 120, 130, 130, 125, 120)

xMusicCurve <- xAudibleCurve[2:(length(xAudibleCurve)-1)]
yMusicCurve <- c(73, 55, 45, 40, 30, 32, 30, 35, 55, 70, 80, 90, 90, 95, 100, 100, 90, 80, 73)

xSpeechCurve <- xAudibleCurve[c(4:9,14:18)]
ySpeechCurve <- c(60, 50, 54, 42, 36, 40, 57, 68, 80, 75, 60)
  
plot_ly() %>% 
    add_trace(x=mat[,1], y=soundIntensity_dB, z=t(m), type="contour", opacity=0.9, # coloring = 'lines'
              line = list(smoothing=0.5, width=3), contours=list(showlabels=TRUE, coloring='heatmap')) %>%
    add_trace(x=xSpeechCurve, y=ySpeechCurve, type="scatter", mode="lines", opacity=0.9, # coloring = 'lines'
              line = list(shape="spline", smoothing=2.5, width=3), hovertemplate = paste('(%{x}, %{y})'),
              name="Threshold of Audibility (bottom) and Pain (top)") %>%
    # add Music curve (parametric spline)
    add_trace(x=xMusicCurve, y=yMusicCurve,type="scatter", mode="lines", opacity=0.9, # coloring = 'lines'
              line = list(shape="spline", smoothing=2.5, width=3), hovertemplate = paste('(%{x}, %{y})'),
              name="Music audible spectrum") %>%
    # add Speech curve domain (parametric spline)
    add_trace(x=xSpeechCurve, y=ySpeechCurve,type="scatter", mode="lines", opacity=0.9, # coloring = 'lines'
              line = list(shape="spline", smoothing=2.5, width=3), hovertemplate = paste('(%{x}, %{y})'),
              name="Speech domain") %>%
    layout(title="Allen & Neely sound-intensity and frequency model",
           xaxis=list(title="Frequency (kHz), log-scale", tickmode = "array", nticks = dim(mat)[1],
                      tickvals = mat[,1], type = "log"), legend = list(x=0.02, y = 0.2), # orientation = 'h'),
           yaxis=list(title="Sound Level (dB)")) %>% #, showlegend = FALSE, showscale = FALSE) %>% 
    hide_colorbar()
# m <- matrix(0, nrow = length(freq_kHz), ncol = length(soundIntensity_dB))
# map_f2i <- function(f) {  # map real frequency (f) to model frequency band (i)
#   i = 1 
#   for (j in 1:(dim(mat)[1]-1)) { # for each of the 11 AN-model frequency bands use appropriate poly-model
#     if (f>=mat[j,1] && f < mat[j+1,1]) { # this is the right model freq-range to get the poly coeffs (mat[,i])
#       i = j
#     } else i = j+1
#   }
#   return(i)
# }
# for (f in 1:(length(freq_kHz)))  {
#   i = map_f2i(f)
#   for (s in 1:length(soundIntensity_dB))  {
#     m[f,s] = mat[i,5] + mat[i,4]*soundIntensity_dB[s] + 
#               mat[i,3]*(soundIntensity_dB[s]^2) + mat[i,2]*(soundIntensity_dB[s]^3)
#   }
# }; image (m)

# Here is an example of a bird chirp sound wave from the package seewave
library(seewave)
## 
## Attaching package: 'seewave'
## The following object is masked from 'package:plotly':
## 
##     export
data(tico)
listen(tico, f=22050)
plot_ly(y=tico@left, type="scatter", mode="lines", name="Tico audio") %>%
    layout(title="Time vs. Amplitude (tico sound)",
           xaxis=list(title="time", tickmode = "array"), 
           yaxis=list(title="amplitude"))
# str(tico)
# See this discussion on how to embed audio player in HTML RMD output:
# https://community.rstudio.com/t/audio-files-in-r-markdown/20874/4 
xTimeSpec <- spectro(tico,f=22050,scale=FALSE, plot=F)$time
yFreqSpec <- spectro(tico,f=22050,scale=FALSE, plot=F)$freq
zAmpSpec <- spectro(tico,f=22050,scale=FALSE, plot=F)$amp
plot_ly(x=xTimeSpec, y=yFreqSpec, z=zAmpSpec, type="contour", opacity=0.9, name = 'Tico Sound',
              line = list(smoothing=0.5, width=0.5), contours=list(showlabels=TRUE, coloring='heatmap')) %>%
    layout(title="Time vs. Frequency Spectrogram (tico sound)",
           xaxis=list(title="time", tickmode = "array"), 
           yaxis=list(title="frequency")) %>% #, showlegend = FALSE, showscale = FALSE) %>% 
    hide_colorbar()
# Because of Heisenberg uncertainty principle, the short-term Fourier transform cannot be localized simultaneously
# in both time and frequency. The temporal and frequency localization dependents on the wl (wave-length) value. 
# High wl values increase the frequency resolution and reduce the temporal resolution and vice-versa, e.g., wl=21
plot_ly(x=yFreqSpec, y=zAmpSpec[,21], 
        type="scatter", mode="lines+markers", name="WL=21 (Tico audio)") %>%
  add_trace(x=yFreqSpec, y=zAmpSpec[,10], 
        type="scatter", mode="lines+markers", name="WL=13 (Tico audio)") %>%
  layout(title="Frequency (Hz) vs. Intensity (dB) (tico sound, at two wavelengths)",
           xaxis=list(title="Frequency (kHz)", tickmode = "array"), 
           yaxis=list(title="Relative amplitude (dB)"),
         legend=list(orientation = 'h'))
# See the Sound/Word-generation RShiny App: https://cogsci.shinyapps.io/soundgen/
# Source code is on Git: https://github.com/tatters/soundgen 

5 Hardy-Weinberg Genotypes Equilibrium

Let’s consider one single nucleotide polymorphism (SNP) in a specific species. Each SNP has two alleles (one of two or more versions of the same gene at the same place on a genetic chromosomal sequence). These alleles could be “called”, or named/labeled, by their nucleotides so that the minor alleles represent the less common versions in the population (allele 1) and the major alleles the more common versions in the population (allele 0).

In genome-wide association studies (GWAS), the minor and major alleles could be labeled differently, as when aggregating data across multiple data sets or populations we need to always double-check the labels for consistency during the preprocessing data harmonization step.

The minor allele frequency (MAF), \(f\), the frequency at which the second most common allele occurs in a given population, is important in heritability. Singletons, MAF variants that occur only once, and more broadly, variants with minor allele frequencies \(\lt 1\%\), drive the vast majority of the observed human mutations.

In most populations, each individual inherits two copies of the genome with each SNP in the sequence having three possible genetic types (genotypes), commonly labeled by the number of copies of allele 1 it contains, i.e. genotype = 0, 1 or 2. Assuming that at a given SNP, the two alleles in one individual are randomly sampled from the entire population, then the the population genotype frequency distribution follows binomial distribution \(\textrm{Bin}(2,f)\).

Genotype Expected frequency from Bin(2,f)
0 \((1-f)^2\)
1 \(2\,f(1-f)\)
2 \(f^2\)

These expected frequencies are called the Hardy-Weinberg equilibrium (HWE) genotype frequencies and define the theoretical equilibrium genotype frequencies given the value of \(f\), the second most common allele frequency. This assume an ideal randomly mating population without selection, migration, or genetic drift.

For most variants in human populations, the observed frequencies approximately follow HWE and deviations from HWE may suggest recent non-random dynamics in the population structure, e.g., mixing of two populations, assortative mating, environmental or other natural selection effects where genotype 1 may be advantageous compared to genotype 2, or one of the genotypes may be lethal become extinct from the population.

In practice, the process of genotyping from the observed intensity on a genotyping chip may have some technical sequencing problems in genotype calling, which can also cause deviations from HWE. This may indicate poor data quality, or alternatively, a variation that is not a bi-allelic, i.e., SNPs may have more than two alleles. Often, variants that do not follow HWE frequencies may be excluded from many GWAS analyses during the quality control (QC) process.

5.0.1 Statistical Test for HWE

To assess the deviation from the HWE, we use statistical inference based on some parametric or non-parametric tests. For instance, we can use a \(df=1\) chi-square distirbution with expected counts derived under the HWE principle with specific allele frequencies.

Suppose the observed genotype counts are \(n_0,n_1,n_2\) for the three genotypes \(0\),\(1\) and \(2\), respectively, and \(N = n_0 + n_1 + n_2\) is the total number of observations. Then, we can estimate the population frequency of allele 1 as \(\hat{f} = \frac{n_1+2\,n_2}{2\,N}\), and the corresponding three expected genotype counts under HWE hypothesis are \(h_0 = N(1-\hat{f})^2\), \(h_1 = 2 N \hat{f}(1-\widehat{f})\), and \(h_2 = N \hat{f}^2\). Then, under HWE assumption,the corresponding \(\chi^2_{df=1}\) test statistic that measures the deviation between the observed counts and the expected counts will be: \[t_{HWE} = \sum_{i=0}^2 \frac{(n_i-h_i)^2}{h_i}\sim \chi^2_{df=1}.\] The value of the test-statistic \(t_{HWE}\) is used for deriving a corresponding \(p\)-value that quantifies the strength of the evidence to reject the null hypotheses (HWE).

Theoretically, the chi-square distribution model is asymptotic and may not perfect for smaller sample sizes or very rare variants. A non-parametric alternative test statistic that does not require this prior assumption is based on permutations distribution, see the R package HardyWeinberg.

5.0.2 Example

Look up the Alzheimer disease SNP rs429358, located in the fourth exon of the ApoE gene on chromosome \(19\). The SNP affects the amino acid at position 130 of the resulting protein. The more common rs429358 allele is (T). When the (C) allele is observed, the same chromosome also harbors the rs7412 (C) allele, and this paired combination is known as an APOE-ε4 allele, which has a strong influence on the risk of developing Alzheimer’s disease.

For the two alleles, \(T\) and \(C\), \(C\) is predicted to be ancestral, the older allele, and \(C\) is the minor allele with average MAF of \(15\%\) across human populations.

The Ensembl Population Genetics shows the allele and genotype frequencies for different human populations. Also, the 1000 Genomes project data for ALL humans shows the minor allele \(C\) frequency \(\frac{|C|}{|T|+|C|}=\frac{754}{4254+754} \approx 15.1\%\) and the observed genotype counts are \(1821\) (TT), \(71\) (CC), and \(612\) (TC) individuals. Note that there are substantial differences between the overall (ALL) and some specific populations, e.g., Northern European descendants, CEU, (T: 0.823 (163); C: 0.177 (35); T|T: 0.667 (66); C|C: 0.020 (2); C|T: 0.313 (31)) and Southern Europeans, TSI, (T: 0.897 (192); C: 0.103 (22); T|T: 0.794 (85); C|T: 0.206 (22); C|T: 0.0 (0)).

We can display the CEU genotype distribution and demonstrate a statistical test for HWE.

# CEU genotype frequencies
genoCEU = c(66, 31, 2)
nCEU = sum(genoCEU) # number of CEU individuals
fCEU = sum(genoCEU*c(0,1,2))/(2*nCEU) #(66*0 + 31*1 + 2*2) / (2*(66+31+2))
fCEU #MAF
## [1] 0.1767677
# Expected genotype frequencies under HWE
hweCEU.prop = c((1-fCEU)^2, 2*fCEU*(1-fCEU), fCEU^2)

# Observed genotype freqs and the HWE
rbind(obs = genoCEU/nCEU, hwe = hweCEU.prop) 
##          [,1]      [,2]       [,3]
## obs 0.6666667 0.3131313 0.02020202
## hwe 0.6777115 0.2910417 0.03124681
# Despite the smaller sample size, this demo uses the HWE Chi-Square test
hweCEU.test = sum((genoCEU - nCEU*hweCEU.prop)^2/(nCEU*hweCEU.prop)) # HWE test statistic
hweCEU.p = pchisq(hweCEU.test, df = 1, lower = FALSE) # p-value corresponding to chi^2 test stats

genotype <- c("GT=0", "GT=1", "GT=2")
plot_ly(x = ~genotype, y = ~genoCEU, type = 'bar', name = 'CEU') %>% 
  layout(title=paste("rs429358 CEU in 1000G Phase3; HWE p=", signif(hweCEU.p,3)),
         yaxis = list(title = 'Proportion'))

Next we will use a larger scale simulation to demonstrate a more realistic GWAS study for the same population of these Northern European descendants, CEU. We will generate the genotype data by first sampling \(n = 10000\) from the genotype frequencies and then sampling independently from the allele frequencies (assuming HWE). As this variant follows HWE, we may not get significant qualitative differences between the two simulations.

# For Rmd reproducibility purposes
set.seed(1234)
n = 10000

# Sample with replacement from genotype frequencies
sampleGeno = sample(c(0,1,2), prob = genoCEU, size = n, replace = T) 

tab = table(sampleGeno) 
# Initialize the carriers of each genotype
countsGeno = rep(0,3)
countsGeno[ 1 + as.numeric( names(tab) )] = as.numeric(tab) #works even if some count is 0

# Sample from HWE frequencies
# sampleHWE = sample(c(0,1,2), prob=c((1-fCEU)^2, 2*fCEU*(1-fCEU), fCEU^2), size=n, replace=T)
# Alternatively, sample n genotypes directly from Bin(2,f) distribution
sampleHWE = rbinom(n, size = 2, p = fCEU)
countsHWE = rep(0,3)
for(i in 0:2) countsHWE[i+1] = sum(sampleHWE == i)

rbind(geno = countsGeno/n, hwe = countsHWE/n)
##        [,1]   [,2]   [,3]
## geno 0.6684 0.3125 0.0191
## hwe  0.6753 0.2942 0.0305
simulation  <- c("Genotype", "HWE")
GT0 <- c((countsGeno/n)[1], (countsHWE/n)[1])
GT1 <- c((countsGeno/n)[2], (countsHWE/n)[2])
GT2 <- c((countsGeno/n)[3], (countsHWE/n)[3])

df <- data.frame(simulation, GT0, GT1, GT2)
df
##   simulation    GT0    GT1    GT2
## 1   Genotype 0.6684 0.3125 0.0191
## 2        HWE 0.6753 0.2942 0.0305
plot_ly(df, x = ~simulation, y = ~GT0, type = 'bar', name = 'GT=0') %>% 
  add_trace(y = ~GT1, name = 'GT=1') %>% 
  add_trace(y = ~GT2, name = 'GT=2') %>% 
  layout(yaxis = list(title = 'Proportion'), barmode = 'stack')

The two simulations appear very similar and we can quantify their statistical differences by assessing the uncertainty of the estimates, e.g., by \(95\%\) confidence intervals for the estimates of each genotype frequency in both of the data sets.

For small counts, a Jeffreys interval provides a Bayesian approach that is better than a freqientist strategy using the standard \(95\%\) confidence interval. For higher counts the two approaches are asymptotically convergent.

# Initialize the matrices
intervalGeno = matrix(NA, ncol = 2, nrow = 3)
intervalHWE = matrix(NA, ncol = 2, nrow = 3)

for(i in 1:3){ # loop over 3 GTs to find intervals
  intervalGeno[i, ] = qbeta(c(0.025,0.975), countsGeno[i]+0.5, n-countsGeno[i]+0.5)
  intervalHWE[i, ]  = qbeta(c(0.025,0.975), countsHWE[i]+0.5, n-countsHWE[i]+0.5)
}

Report the observed genotype frequencies (column=1) and its \(95\%\) Jeffreys interval (columns 2 and 3) for both data sets. Then compare whether the differences between the estimates using the computed uncertainty about the intervals. The point estimates are within other data set’s 95% Jeffreys credible interval and the frequency differences between the two data sets are not substantive. Again, we can use the standard two-sample \(\chi^2\) test to quantify the frequency difference as a probability values.

cbind(genoEst = countsGeno/n, genoEst_LB=intervalGeno[,1], genoEst_UB=intervalGeno[,2],
      hweEst = countsHWE/n, hweEst_LB=intervalHWE[,1], hweEst_UB=intervalHWE[,2]) 
##      genoEst genoEst_LB genoEst_UB hweEst hweEst_LB hweEst_UB
## [1,]  0.6684 0.65912503 0.67757751 0.6753 0.6660725 0.6844260
## [2,]  0.3125 0.30347059 0.32163793 0.2942 0.2853293 0.3031898
## [3,]  0.0191 0.01655458 0.02192369 0.0305 0.0272643 0.0340074
chisq.test(rbind(countsGeno, countsHWE))
## 
##  Pearson's Chi-squared test
## 
## data:  rbind(countsGeno, countsHWE)
## X-squared = 32.076, df = 2, p-value = 1.084e-07

6 Levers

7 Motion

8 Stress, Strain & Sheer

9 Strings and Dynamic Compression & Stretching

10 Energy & Fracturing

11 W Bone fracture

11.1 Impulses

11.2 Falls

12 Work

13 Compressibility and Viscosity

14 Human Circulatory System

SOCR Resource Visitor number Web Analytics SOCR Email