Methods and systems for estimating blood analyte conditions

Abstract

Methods and systems for estimating blood analyte conditions. In some methods, signal data may be received from a non-invasive blood monitor. The signal data may then be processed to remove a mean offset component from the signal data. The processed signal data may be transformed to a frequency domain. A DC offset component may then be extracted from the transformed signal data, after which the remaining data, which may comprise data from which the extracted data was taken, may be used to predict a blood analyte condition.

Description

METHODS AND SYSTEMS FOR ESTIMATING BLOOD ANALYTE CONDITIONS CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit under 35 U.S.C. § 119(e) of U.S. Provisional Patent Application No. 63 / 601,494, which was filed November 21, 2023, and titled “Integrated Path Approximation: A Novel Approach to Modeling Macroscopic Photon Absorption in Highly Scattering Media,” which is hereby incorporated herein by reference in its entirety. SUMMARY

[0002] Noninvasive blood monitoring is a decades-old research problem that has proved to be frustratingly elusive. Data generated from a state-of-the-art infrared spectrometer, combined with new techniques to model physiological interactions, may reveal previously inaccessible signals. Disclosed herein are methods, which may be incorporated into various systems, for deriving distributional estimates for photon scattering and absorption behavior, which lead to simple closed-form models that accurately model photon scattering. These analytic models improve reliable physiological signal extraction and continuous noninvasive measurement of biomarkers of interest.

[0003] To reliably determine small signals from noisy data, reasonable models of the underlying signal-generating process must be developed. This has proved to be difficult in the case of, for example, glucose measurement, notably because the light-particle interactions present are difficult to model in a closed-form way. Monte Carlo methods exist, but such Monte Carlo methods would be intractable for real-time measurements on embedded hardware.

[0004] Disclosed herein are closed-form models that reasonably approximate light- particle interactions in the body, and may therefore facilitate data assimilation efforts, along with estimation and / or prediction of various conditions, such as blood analyte conditions. Methods for extracting DC offset data from a data set are also disclosed herein. In some embodiments and implementations, such offset data may be extracted from a data set comprising a pulsatile signal and a non-pulsatile blood component. The DC offset component in such cases may comprise data derived from tissue, which may 1 Docket No. JTS0005be extraneous for purposes of analyzing analytes in the blood. Such extracted data may therefore, in some cases, be discarded. Extracting such data may therefore be important in enhancing the accuracy of data processing techniques intended to estimate, predict, and / or correlate data with blood analyte conditions.

[0005] Techniques for removing and / or extracting certain components from a pulsatile signal used for estimating certain conditions, such as blood analyte conditions, are also disclosed. For example, because the DC offset component of a non-invasive electromagnetic radiation based sensor system typically contains data that may not be useful, and may make it more difficult to obtain useful information from the data, techniques disclosed herein may be used to extract the DC offset component, which in the case of a non-invasive blood monitoring system may leave only data associated with the pulsatile signal (e.g., the heartbeat being transduced through a blood vessel wall), along with data corresponding to the non-pulsatile blood component of the signal (i.e., electromagnetic waves reflected from the blood itself).

[0006] In an example of a method for processing signals from a non-invasive blood analyte tracking system, the method may comprise receiving signal data from a non- invasive blood monitor. The signal data may comprise a frequency component corresponding to a heart rate pulse of a user of the non-invasive blood monitor; a non- pulsatile blood component corresponding to light reflected from within a blood vessel; and a tissue-dependent DC offset component. The signal data may be processed to extract one or more components, such as the tissue-dependent DC offset component. The processed signal data may then be evaluated to provide an estimate of a blood analyte condition using only one or more of the remaining components, such as the frequency component and the non-pulsative blood component, from the signal data.

[0007] In some implementations, the blood analyte condition may comprise a concentration of the blood analyte.

[0008] In some implementations, the blood analyte condition may comprise a trend over time associated with the blood analyte.

[0009] In some implementations, the blood analyte may comprise glucose. Alternatively, or additionally, the blood analyte may comprise oxygen.

[0010] In some implementations, the step of processing the signal data may comprise 2 Docket No. JTS0005use of a probability density model, such as a Levy Distribution.

[0011] In an example of a method for extracting a DC offset component from signal data using a non-invasive blood monitor according to some implementations, the method may comprise receiving signal data from a non-invasive blood monitor and processing the signal data to remove a mean offset component from the signal data. The processed signal data may be transformed to a frequency domain. A DC offset component may then be extracted from the transformed signal data. The remaining data, such as pulsatile signal data and / or non-pulsatile signal data from the blood, may then be used to predict a blood analyte condition.

[0012] In some implementations, the blood analyte condition may comprise a concentration of the analyte and / or a trend over time associated with the analyte.

[0013] The blood analyte may comprise any desired constituent of blood, such as glucose or oxygen, for example.

[0014] Some implementations may further comprise adjusting signal data to reduce scale discrepancies.

[0015] Some implementations may further comprise processing the signal data to compute discrepancy measure.

[0016] Some implementations may further comprise adjusting the offset to reduce discrepancy.

[0017] In an example of a method for estimating a blood analyte condition using data from a non-invasive blood monitor according to some implementations, the method may comprise receiving signal data from a non-invasive blood monitor. The signal data may comprise various separable constituents. For example, in some cases, the signal data may comprise a frequency component corresponding to a pulsatile signal, such as a heart rate pulse, of a user of the non-invasive blood monitor; a non-pulsatile blood component corresponding to light reflected from within a blood vessel; and a tissue-dependent DC offset component. One or more of these components may then be extracted. For example, some implementations may comprise extracting the non-pulsatile blood component and / or the DC offset component from the signal data. The extracted component(s) may then be used to estimate a blood analyte condition, such as a blood glucose level or blood oxygen level, for example. 3 Docket No. JTS0005

[0018] Some implementations may further comprise extracting the frequency component. In some such cases, the frequency component may be used with the non- pulsatile blood component to estimate the blood analyte condition.

[0019] The blood analyte condition may comprise, for example, a trend over time in a concentration of the blood analyte and / or a concentration of the blood analyte.

[0020] The features, structures, steps, or characteristics disclosed herein in connection with one embodiment may be combined in any suitable manner in one or more alternative embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] The written disclosure herein describes illustrative embodiments that are non- limiting and non-exhaustive. Reference is made to certain of such illustrative embodiments that are depicted in the figures, in which:

[0022] FIG.1 is a diagram illustrating a, arrangement for a standard Beer-Lambert law experiment;

[0023] FIG.2 is a diagram illustrating a single scattering event in which ϕ is the change in rotational angle from the old path, and θ is the change in azimuth angle from the old path;

[0024] FIGS. 3A and 3B are graphs comparing the Von Mises Distribution with the wrapped normal distribution;

[0025] FIGS.4A and 4B are graphs comparing the density and cumulative density of the Von Mises vs. Uniform Distributions;

[0026] FIGS.5A and 5B are graphs illustrating True CDF vs. Arctan Estimate of HG function, along with samples generated using inversion sampling;

[0027] FIG. 6A is a graph illustrating the Hoeffding dependence across scattering probabilities;

[0028] FIG.6B is a graph illustrating the density of samples taken;

[0029] FIG. 7 is a graph illustrating an empirical CDF of samples generated from Monte-Carlo simulation of photon update equations according to some implementations compared to the modified Levy distribution;

[0030] FIG.8 is a diagram illustrating a framework for a model involving photon paths 4 Docket No. JTS0005from the surface of a user’s skin to an artery underneath the skin;

[0031] FIG. 9A is a graph plotting experimental data using the loss function and plotting the loss for various values of x and y corresponding to various DC offset values;

[0032] FIG.9B is another graph plotting synthetic data generated using the techniques of certain implementations described herein;

[0033] FIG.10A is a graph schematically illustrating the constituent components of a signal having a pulsatile portion, a non-pulsatile blood component, and a DC offset / tissue component;

[0034] FIG. 10B is a diagram illustrating the various electromagnetic waves that correspond to the signal components illustrated in FIG.10A;

[0035] FIG.11 is a flowchart illustrating a method for extracting a DC offset component of a signal according to some implementations;

[0036] FIG. 12 is a block diagram illustrating an example of a system for extracting one or more components from a pulsatile signal; and

[0037] FIG.13 is a flowchart illustrating an example of a method for estimating a blood analyte condition according to some implementations. DETAILED DESCRIPTION

[0038] It will be readily understood that the components of the present disclosure, as generally described and illustrated in the drawings herein, could be arranged and designed in a wide variety of different configurations. Thus, the following more detailed description of the embodiments of the apparatus is not intended to limit the scope of the disclosure, but is merely representative of possible embodiments of the disclosure. In some cases, well-known structures, materials, or operations are not shown or described in detail.

[0039] As used herein, the term “substantially” refers to the complete or nearly complete extent or degree of an action, characteristic, property, state, structure, item, or result to function as indicated. For example, an object that is “substantially” cylindrical or “substantially” perpendicular would mean that the object / feature is either cylindrical / perpendicular or nearly cylindrical / perpendicular so as to result in the same or nearly the same function. The exact allowable degree of deviation provided by this term may depend on the specific context. The use of “substantially” is equally applicable when 5 Docket No. JTS0005used in a negative connotation to refer to the complete or near complete lack of an action, characteristic, property, state, structure, item, or result. For example, structure which is “substantially free of” a bottom would either completely lack a bottom or so nearly completely lack a bottom that the effect would be effectively the same as if it completely lacked a bottom.

[0040] Similarly, as used herein, the term “about” is used to provide flexibility to a numerical range endpoint by providing that a given value may be “a little above” or “a little below” the endpoint while still accomplishing the function associated with the range.

[0041] The embodiments of the disclosure may be best understood by reference to the drawings, wherein like parts may be designated by like numerals. It will be readily understood that the components of the disclosed embodiments, as generally described and illustrated in the figures herein, could be arranged and designed in a wide variety of different configurations. Thus, the following detailed description of the embodiments of the apparatus and methods of the disclosure is not intended to limit the scope of the disclosure, as claimed, but is merely representative of possible embodiments of the disclosure. In addition, the steps of a method do not necessarily need to be executed in any specific order, or even sequentially, nor need the steps be executed only once, unless otherwise specified. Additional details regarding certain preferred embodiments and implementations will now be described in greater detail with reference to the accompanying drawings.

[0042] A classical way to understand light-particle interactions is the Beer-Lambert law, which states that the macroscopic behavior of light when a narrow-beam light source propagates through a diffusive media is governed by the following equation: ^^ = ^^0^^−^^^^^^

[0043] FIG. 1 illustratesfor the standard Beer-Lambert law experiment and equation, in which ^^ is the intensity of detected / transmitted light, ^^0is the initial intensity of the light, ^^ is the molar absorption coefficient (material dependent), ^^ is the optical path length, and ^^ is the concentration of the absorptive species in the media. Notably, this equation only applies when the path length is constant for all detected paths (or can be approximated to be constant) and the media has a uniform distribution 6 Docket No. JTS0005of absorbing media.

[0044] In the case of a body-spectrometer system, however, such approximations do not reflect the physical reality of the system. Other work related to physiological signal extraction includes the differential Beer-Lambert law, useful in continuous-wave near- infrared tissue spectroscopy (cwNIRS) for measuring changes in tissue chromophore concentrations, particularly oxy- and deoxyhemoglobin. The differential Beer-Lambert law is expressed as: ^^^^ = ln^^det,1^^= ^^^^^^^^

[0045] where ^^^^ represents attenuation, ^^det,1and ^^det,2are thedetected light intensities in two different tissue states, ^^ is the mean path length of detected photons, and ^^^^^^is the change in the absorption coefficient.

[0046] This law posits that the change in light attenuation is directly proportional to changes in chromophore concentrations, allowing for the quantification of these changes based on detected light intensities. However, the differential Beer-Lambert law relies on two critical assumptions: (1) homogeneous absorption changes within the tissue, and (2) constant scattering loss. Deviations from these assumptions, such as inhomogeneous absorption changes and variations in scattering, can lead to underestimations in concentration changes (partial volume effect) and crosstalk between chromophores..

[0047] This disclosure provides a new model that is more expressive than the differential model, as it models not only the change in absorption but also absolute absorption values, thereby providing more information about the observed physiological system. Moreover, the model derived herein is more representative of the true physics of photon scattering in turbid media, which may therefore provide a suitable path for improving signal processing in blood and other body fluids.

[0048] Instead of considering all photons simultaneously, it may be useful to study the path of a single photon packet at a time. We can consider a small part of the path of this photon and may be able to make several assumptions. For example, we may assume that this section of path is small enough that the scattering media can be assumed to be uniformly distributed along the section. We may also be able to assume that this section of the path is large enough that we still observe the macroscopic effects of the Beer- 7 Docket No. JTS0005Lambert law. Finally, we may assume that refraction due to Snell’s law is negligible given the multitude of other scattering events.

[0049] Given these assumptions, we can model such a section of the path as: ^^^^+1 = ^^^^^^−^^^^^^^^^^

[0050] where ^^^^and ^^^^depend on time, and ^^ is fixed. Since this is a recursive equation, we can write: ^^^^+1 = ^^0^^∑^^^^=0 − ^^^^^^^^^^

[0051] Since each path is a absorbing barrier (i.e., absorbed ifthe particle leaves the body or is , we may consider the final time ^^^^^^. Sincesummation is linear, and by the property of ^‾^ the algebraic mean ∑^^ ^^ ^^^^ = ^^^‾^, we canwrite the following: ∑^^^^ − ^^ ^^^^ − ∑^^ ^^^^^^ = ^^0^^^^=0^^ ^^ = ^^0^^ ^^ ^^ ^^^^ ^^^^ = ^^0^^−^^^^^‾^

[0052] where A isleaving the body orbeing detected, and ^‾^ is the average absorption along the path. Note that ^‾^^^ is combinedwith ^‾^^^ since they are indistinguishable in the sum. Thus, we have reduced the absorptionalong a path down to a Beer-Lambert like expression, but this only applies to a single photon path. In order to calculate the resultant detected light at a particular wavelength, we have the expression: ^^ ^^ ^^0 ^^−^^^^^^^‾^^^

[0053] where ^^ is the bulksum over detected paths), ^^0is total initial intensity of light, (thus^^0^^ is the intensity of light per packet,), ^^^^is the numberof steps the photon packet took before being detected, and ^‾^^^ is the average absorptionof a single path.

[0054] It is worth noting the similarity of this expression to the Feynman-Kac formula, exchanging integrals for summation, as necessary. This is an interesting connection that bears more investigation. Now, if we wish to simplify this model further, we must understand ^^^^and ^‾^^^. 8 Docket No. JTS0005

[0055] We begin with ^‾^^^, as it is the easier of the two. We can consider each ^^^^as a sample from a categorical distribution since for each step along the path, we end up inone of ^^ different body tissue compartments, with some probability ^^1... ^^^^. Note that wemay consider this categorical distribution to be fixed across one path, since the summation in the recursive approximation is linear, so we can rearrange tissue compartments as we like.

[0056] Proposition 1: ^^ −^^ ∑^ ^^^^ = ^^0^^ ^^^^^^^ ^^^^ = ^^ −^^ ∑^^ 0^^^^^̃^^^^^^^where the ^^^̃^are i.i.d. categorical random variables defined only

[0057] Proof: ^^^‾^ = ∑^^ ^^=1 ^^^^and we have that given ^^ different possible values of ^^, thensince ^^[^^] = ∑^^ ^^=1 ^^^^^^^^, we can define a categorical distribution to make it work.

[0058] we have ^^^‾^ = ∑^^ ^^ ^^^^, and since all we observe is the bulk=1absorption across the path, we can reorder the ^^^^as we like, thus we can take them as i.i.d. Thus, each ^^^^is a sample from some categorical distribution, with each ^^^^independent from the others. Thus the ^^^^are i.i.d. and so by the Central Limit Theorem(CLT), the ^‾^ have a normalNote this distribution is not necessarily thedistribution across different paths; in fact, we expect it to be different.

[0059] To understand the ^^^^, we must embark on a journey of probability theory. Wecan first re-frame our definition of ^^^^. We have that ^^^^^^ = ^^ for some total path length L. Ifwe instead consider small time steps over our photon path and assume the velocity of light is constant across the path (a reasonable approximation since the light is mainlypropagating through water, no matter the tissue compartment), then we can write ^^ = ^^^^^^,where ^^^^is now the time the photon has spent in the body, and ^^ is the speed of light in the body. Given this assumption, we will henceforth refer to a path’s length in terms of its time since the two are equivalent given constant velocity.

[0060] In preferred embodiments, an important feature of the model is the fact that the ^^^^are distributed according to a modified Levy Distribution, which is the distribution of hitting times of Brownian motion. This section is dedicated to the proof of that fact. Thedensity function for the hitting time of the ^^^^ will henceforth be denoted ℒ( ⋅   | ^^, ^^), definedas: 9 Docket No. JTS0005ℒ(^^|^^, ^^) =^^−3 / 2 ^^2 / (2^^^^2)√ ^^^^2^^ ^^ , 2

[0061] where ^^ is the to the variance of therandom walk, or rather the , the distance of the absorbing (detection) barrier from the starting point of the path. This distributional information will allow us to make simplifications that respect the underlying probability laws, which provides simple, closed-form approximations that are accurate and computationally efficient. Much of the work done here is to develop and justify these simplifications.

[0062] To prove that the density of the ^^^^is given by the equation above, we can simplify the three-dimensional case to the one-dimensional case. We note that to understand when paths leave the skin, we need only consider the movement in the z- direction, where we consider the x and y directions to be coplanar with the surface of the skin. We also consider the skin to be a flat Cartesian plane, ignoring the curvature of the surface of the skin. For most purposes, this is a reasonable approximation since any photon detector is localized to a small area of the wrist, where curvature is typically negligible. Moreover, any photons that travel a significant distance from the incidence point of the photon source will have been mostly absorbed, so any contributions from these photons will be negligible.

[0063] Thus, to understand the distribution of lengths of paths, we may simply consider the length of time they spend in the tissue, which will be determined by the time they leave the surface of the skin. Then, the expression we must consider is only the z component of the random walk of the photons.

[0064] We can express the position of a photon using standard spherical coordinates, with ^^ = ^^sin(^^)cos(^^), ^^ = ^^sin(^^)sin(^^), ^^ = ^^cos(^^). Then, scatteringevents will change the direction according to ^^ and ^^. We may then consider such a scattering event as illustrated in FIG.2, where ^^ is the change in inclination from the old trajectory, and ^^ is the change in rotation angle from the old trajectory.

[0065] A single scattering event illustrated in FIG. 2, in which ϕ is the change in rotational angle from the old path, and θ is the change in azimuth angle from the old path. The old path is illustrated with the dotted line, and the dashed lines illustrate the change in angle in the x,y and y,z planes, respectively.Docket No. JTS0005

[0066] Then, in order to consider the absolute position, we must add the update angles, since if a trajectory is moving at some angle ^^1with respect to the skin, (or some fixed reference point) and ^^2is the change in angle from the original trajectory, then thetotal angle, relative to the skin will be ^^1 + ^^2, modulo 2^^. Then, the absolute position ofthe photon at a given time step will be given by the following update equations: ^^ ^^ ^^^^+1 = ^^^^ + ^^ ⋅ sin (∑ ^^^^ ) cos (∑ ^^^^) )

[0067] Then, we can z as a random walk, and we consider updates to the z position as drawn from some distribution. Photon scatteringis often modeled using a Heyney-Greenstein distribution, denoted ^^^^(⋅ |^^, ^^), whichaccounts for the anisotropy of scattering, parameterized by ^^, and the contribution of Rayleigh scattering, parameterized by ^^. This density function gives the probability of scattering at a certain angle ^^ away from the current trajectory of the photon. Then, for a typical photon scattering event, we will have the update angle drawn from the modified Henyney-Greenstein distribution, a standard photon scattering model in turbid media. Here we use a slightly modified distribution called the Draine distribution, which accounts for the effects of Rayleigh scattering, parameterized by ^^. We then define the density ^^^^: 31 − ^^2 1 + ^^cos2(^^)^^ =

[0068] To calculate the updates to both ^^ and ^^ for a scattering event, we may assume that once the Heyney- Greenstein angle ^^ is determined, the scattering direction will be uniformly distributed in the "cone" around the original trajectory, with vertex angle ^^. Thus, we will have that, for 11 Docket No. JTS0005a given timestep, the updates will be given by ^^^^ = ^^^^sin(^^^^) and ^^^^ = ^^^^cos(^^^^), where^^^^is uniformly sampled over[0,2^^].

[0069] Now, we note that scattering events will occur with some likelihood ^^^^, such that over a fixed distance, a photon will have some probability ^^^^of scattering or not. Thus, if we wish to update the absolute orientation of the photon, we will update the positions by multiplying ^^^^by a random variable sampled from a Bernoulli distribution with likelihood ^^^^.

[0070] Thus, the absolute z position of the particle at time index ^^ + 1 will be given by:^^ ^^ ^^^^+1 = ^^ ⋅ cos ^^^^sin

[0071] Where ^^^^   ∼ ^^^^   ∼  ^^^^(^^, ^^).

[0072] The independence of scattering events can be modeled as follows: ^^ ^^ =^^ ⋅ ^^^^from^^ ⋅^^ ⋅independent, andwalk with the strong Markov property

[0073] Thus, to model this position as a random walk with the strong Markov property, we wish to find some step size ^^ so that we can express the updates to the ^^ position as independent, identically distributed (i.i.d.) random variables. This is justified since the ^^^^sin(^^^^)^^^^are i.i.d., and thus their sum will converge to a normal distribution by CLT. Then in order to get independence, we must pick ^^ large enough so that after ^^ steps, 12 Docket No. JTS0005the positions become effectively independent, allowing us to encapsulate the ^^ updatesas a single random variable, which will be independent of the ^^ + 1 to 2^^ positionupdates. This independence fundamentally comes from the aliasing of the cos(^^)functionover its domain. Since cos(^^) = cos(^^ + 2^^), then the aliasing of the domain "wraps" thedistribution around the unit circle, a notion formalized with the wrapped normal distribution. The wrapped normal distribution has density function given by: ∞ 1−(^^ − ^^ + 2^^^^)2^^^^^^=∑ exp( )

[0074] where ^^, ^^ normal distribution. As^^ approaches infinity, a on [−^^, ^^]. Thensince the sum of the ^^^^sin(^^^^)^^^^will have a normal distribution with strictly increasingvariance, then the cos(^^) will result in samp ∑∞les from cos( ^^ ^^^^sin(^^^^)^^^^)approaching cos(Uniform(−^^, ^^)).

[0075] Then to getif we pick ^^ large enough, we will have∑^^ ^^=0 ^^^^sin(^^^^)^^^^normally distributed with some large variance. Then the wrappeddistribution will be near to Uniform(−^^, ^^), as will ∑2^^ ^^=^^ ^^^^sin(^^^^)^^^^. Then we canapproximate these two sums as being independentwill have for any ^^ ∈ [−1,1]:^^ 2^^ ^^ ^^

[0076]

[0077] Lemma 1: The sum of two uniform random variables supported on the unit circle will be uniformly distributed.

[0078] Then, if we label: ^^

[0079] then we will have13 Docket No. JTS0005^^(cos(^^1 + ^^2) ≤ ^^) ≈ ^^(cos(^^2) ≤ ^^),

[0080] which implies ^^(cos(^^1) ≤ ^^, cos(^^1 + ^^2) ≤ ^^) ≈ ^^(cos(^^1) ≤ ^^, cos(^^2) ≤ ^^)≈ ^^(cos(^^1) ≤ ^^)^^(cos(^^2) ≤ ^^).

[0081] The second line follows from the fact that shifts thebase point of the distribution of ^^2. Thus, (^^) over itsdomain has ^^2 mod 2^^ distributed nearly uniformly, so by the aforementioned Lemma, we have that their sum is also uniformly distributed.

[0082] Thus, knowing the value of cos(^^1) does not influence the value of cos(^^1 + ^^2),due to the density of ^^2 mod 2^^ being uniformly constant. Therefore, we have thatcos(^^1 + ^^2) and cos(^^2) areand so we have that for some ^^^^ ^^ ^^ ⋅ ^^^^

[0083] is independent from2^^ ^^ ^^ ⋅ ^^^^

[0084] Thus, we can write:^^^^(^^+1) = ^^^^⋅^^ + ^^^^

[0085] Where the ^^^^are

[0086] Atest for independence of cos(^^1) and cos(^^1 + ^^2), with ^^ given bythe lower bound, we derive below. This test demonstrates empirical support for the validityof the assumptions of cos(^^1), cos(^^2) for sufficiently large M.

[0087] In order to continue our analysis on the ^^^^, we will want ^^ to be as small as possible while still retaining independence. Thus, we derive a lower bound for ^^.

[0088] We need to work with the wrapped normal distribution to extract this lower bound. However, the wrapped normal distribution is difficult to work with, so we instead work with a good approximation to the wrapped normal distribution, the von Mises distribution. The wrapped normal distribution describes a freely diffusing angle ^^ on a circle with an unwrapped variance that grows linearly in time. On the other hand, the von 14 Docket No. JTS0005Mises distribution is the stationary distribution of a drift and diffusion process on the circle in a harmonic potential, i.e., with a preferred orientation.

[0089] The density function for the von Mises distribution is given by: exp( ( )( ^^cos ^^ − ^^ )^^ ^^|^^, ^^) =

[0090] where ^^0(^^)is the first kind of order 0. Here, ifa von Mises distribution is a normal distribution, we have^^ =1 ^^2, where ^^2is the variance parameter of the wrapped normal distribution, and ^^ is the as in the wrapped normal distribution.

[0091] In FIGS. 3A and 3B, we show plots of the von Mises and wrapped normal distributions, respectively, to illustrate their similarity. It may be preferable to use the von Mises distribution because it is typically simpler to work with.

[0092] Then, in order to isolate a lower bound, we first set a target value ^^ for the vonMises distribution. if ^^ =1 ^^2, then the corresponding von Mises distribution is very nearlyuniform. FIGS. 4A and 4B provide a comparison of the von Mises distribution with ^^ =1 ^^2and uniform (−^^, ^^). The Wasserstein distance between the two is 0.066.

[0093] This value for ^^ corresponds to a variance ^^ = ^^ for the corresponding normaldistribution. This value of ^^ heuristically comes from the fact that the derivative of the normal pdf will have a maximum when the second derivative is zero and a negative third derivative. We wish to characterize the maximal rate of change of the normal density function since we want to approximate a constant distribution once we wrap the normal distribution. Thus, by characterizing the maximal change, we can bound the uniformity of the wrapped distribution by the maximal rate of change of the unwrapped distribution.

[0094] The derivative of the normal pdf is given by: ^^2^^ ^^^^−2^^2

[0095] Where we have setof ^^ is mean centered. Then,the derivative will be maximized when the second derivative: 15 Docket No. JTS0005^^2 −^^2 2^^2(^^2 − ^^2)^^^^ (^^|^ )^^ ^^2^^ ^, ^^ =√2^^^^5

[0096] is zero, or^^3^^2^^2 (−3^^2^^ + ^^3)^^^^( ^^) =

[0097] is negative at ^^

[0098] Then, sincecos wrap distribution onto [−^^, ^^], wehave if ^^ = ^^, then the maximal rate of change of the underlying distribution will bewrapped to the endpoints of the support of the wrapped distribution, and we will have some of this change offset by the wrapping at the point of maximal change. Then thecorresponding von Mises distribution with ^^ =1 ^^2is a reasonable approximation of the uniform distribution, and thus our results relying on the invariance of thefinal distribution of cos(^^) to the aliasing split location go through.

[0099] Then we wish for ∑^^ ^^ ^^^^sin(^^^^)^^^^to be normally distributed with variance ^^2. Thus, we analyze the^^^^in order to find a lower bound on ^^. Wehave that ^^^^ is Bernoulli, which is discrete but can be represented by ^^^^^^(^^ − 1) +(1 − ^^^^)^^(^^), where ^^(^^) is the Dirac distribution. Then ^^^^^^^^ will be distributed accordingto: ^^^^^^^^ ∼ ^^^^^^^^(^^|^^, ^^) + (1 − ^^^^)^^(^^),

[0100] since if ^^^^ = 1,0, then ^^^^^^^^ = 0. Then, sinceeach of these events has probability ^^^^ and 1 − ^^^^, respectively, we simply weight eachof the two cases by ^^^^, and (1 − ^^^^). Then the final distribution will be the weighted sumdue to the independence of ^^^^and ^^^^, so we have by the law of total probability: ^^(^^^^^^^^ ≤ ^^) = ^^(^^^^^^^^ ≤ ^^|^^^^ = 1)^^(^^^^ = 1) + ^^(^^^^^^^^ ≤ ^^|^^^^ = 0)^^(^^^^ = 0)of ^^^^^^^^will be given by: 16 Docket No. JTS0005^^ ^^^^^^(^^^^^^^^) = ∫ ^^2(^^^^^^^^(^^|^^, ^^) + (1 − ^^^^)^^(^^)) =−^^

[0102] since the

[0103] Then we the distribution for sin(^^). In order to use this transform, we need a bijective mapping so the inverse function is well-defined. Then, since sin(^^) is an odd function, we can consider ^^^^distributeduniformly on [0, ^^], and then the distribution for the transform of a uniform random variableon [−^^, ^^] will be an appropriately reflected and scaled version of the distribution forthat sin(^^) is not bijective on [0, ^^] but sin(^^) is symmetric on this domainabout ^^ / 2, and so we can simply consider sin(^̂^^^) for ^̂^^^ distributed uniformly on [0, ^^ / 2].Then by the probability integral transform, ( )Uniform 0, ^^ / 2 :^^(sin(^̂^^^) ≤ ^^) = ^^ (^̂^^^ ≤ ^^ ) =2arcsin(^^)( )^^

[0104] For ^^ ∈ [0,1]. Then we have by theof calculus:2arcsin(^^) ^^ ^^ 2 ( ) ^^ ^^(sin(^̂^^^) ≤ ^^) = = ∫arcsin ^^ ^^^^ = ∫2 ^^^^ ^^ ^^^^ ^^^^

[0105] reflect thedistribution across the y axis and scale appropriately so the distribution integrates to unity. 2 However,^^√1−^^2is symmetric about ^^ = 0, thus our density for sin(^^^^) is given by:^^ (^^) =1

[0106] for ^^ ∈ [−1,1]

[0107] Which has mean 0, (by symmetry), and variance given by: 12^^^^^^ (sin =^^ 1 = 2

[0108] Then we have using the^^(^^2^^2) + ^^(^^^^)2, that thevariance of ^^^^^^(sin(^^^^)^^^^^^^^) is given by:17 Docket No. JTS0005^^^^^^(sin(^^^^)^^^^^^^^) =2 2 ^^^^^^ (sin(^^^^)) ^^^^^^(^^^^^^^^) + ^^^^^^ (sin(^^^^)) ^^(^^^^^^^^) + ^^ (sin(^^^^)) ^^^^^^(^^^^^^^^) =

[0109] ofvariance computation, if we assume ^^^^are normally distributed, we haveby the property of normal distributions ^^(0, 1) + ^^(0, ^^22) = ^^(0, ^^21 + ^^22), that: ^^ ∑^^^^ sin(^^^^)^^^^^^^^^^ ∼ ^^ (0, ^^^^^^(^^^^(^^|^^, ^^)))

[0110] Thus, if ^^2, then we set:2^^2^^ =)

[0111] In order to samplewe must sample the ^^^^from the ^^^^ distribution. However, we wish to perform inversion sampling (so that samples can be generated quickly for other, more complex Monte-Carlo experiments), we must find aninverse to the CDF of ^^^^(^^|^^, ^^). However, this CDF does not have an easily computableinverse thus, we use a parametric function of the form ^^ + ^^arctan(^^^^), and then giventhe CDF of an ^^^^ distribution (which we can compute analytically), we can then fit a least squares estimate of the CDF, and then invert the parametric equation to find ^^−1for^^^^(^^|^^, ^^). We illustrate the validity of this approach in FIGS. 5A and 5B for ^^ = 0.6, ^^ =0.1, which are estimated parameters of the anisotropy of scattering in tissue by Tseng et al. (S. H. Tseng, A. Grant, and A. J. Durkin, In vivo determination of skin near-infrared optical properties using diffuse optical spectroscopy, Journal of Biomedical Optics, 13(1):014016, Jan-Feb 2008.). These figures illustrate the true CDF vs. approximate CDF of HG function, along with samples generated using inversion sampling.

[0112] In some preferred implementations, the next step is to perform Hoeffding tests of independence across multiple values of ^^^^. Here we have set ^^1=∑^^^^ ^^ ^^=1 cos (∑ ^^=1 ^^^^ sin(^^^^)^^^^) and ^^2 = ∑2^^ ^^=^^ cos (∑ ^^=1 ^^^^sin(^^^^)^^^^).Docket No. JTS0005

[0113] Then, for each value of ^^^^, we can sample 2000 pairs ^^1, ^^2 and compute theHoeffding coefficient for the two sample sets. Two independent variables will have a Hoeffding coefficient of zero, while negatively correlated variables will have a coefficient of −1 / 2, and positively correlated variables will have a coefficient of 1. This test will give empirical support for the independence of the ^^^^. We may also analyze the sum of the ^^^^for normality. Here we have sampled 5000 values of ^^1^^^^, where ^^1^^^^ = ∑^^ ^^=1 ^^^^, and ^^ is chosen so that ^^^^ corresponds to 1 nanosecond of photon travel. Thus, given ^^ according to (30), we have ^^ = 2.25 × 102 / ^^, where the numerator corresponds to thenumber of millimeters light travels in water per nanosecond. Thus, ^^1^^^^will be the ^^ position of a photon after 1 nanosecond of time, given the photon travels in water.

[0114] FIG.6A is a plot of Hoeffding coefficients for independence and FIG.6B is a histogram of the ^^1^^^^ for ^^^^ = 1. These figures illustrate the independence of the ^^^^measured using Hoeffding coefficients, along with the density of samples from ^^∑ ^^ ^^=1 ^^^^. We can see that the ^^^^are likely to be independent, and it appears that T is large

[0115] We can then test for normality of the ^^^^using the Shapiro-Wilks test. This test gives empirical support for the convergencethe ^^^^to the normal distribution by CLT, which might not be a good assumption if the number ^^^^required to reach normality is very large (if M is very large, then T will be small, and so assuming the ^^^^are normally distributed may not be justified.). The results below indicate that ^^ is indeed large enough that the ^^^^may be assumed to be normal, at least for values of ^^^^in the range[0.1,1], as indicated by Tseng et al. The Shapiro-Wilks test will producep-value depending on the normality of the sample. If the p-value is large, we fail to reject the null hypothesis that the samples are normally distributed, which indicates that the samples might be normally distributed.

[0116] Below is a table of the median Shapiro-Wilks score (denoted SW) across ten tests for various ^^^^in [0.1,1]. This table includes values for ^^^^and the median Shapiro- Wilks test p-Here, we see large p-values across all of the ^^^^, thus failing to indicate that indicating that the^^1^^^^are not drawn from a distribution that is not likely to be statistically significantly different from a normal distribution. 19 Docket No. JTS0005^^^^0.1 0.14 0.2 0.29 0.42 0.6 0.86 1.0 SW 0.425 0.684 0.595 0.373 0.122 0.201 0.131 0.188

[0117] Given these results, we see no reason we cannot assume that the ^^^^are i.i.d, and the ^^^^are normally distributed. Therefore, our lower bound M allows to continueour analysis, and we can write: ^^^^+1 = ^^^^ + ^^^^

[0118] Where the ^^^^are independent and are given by: ^^ ^^^^^^ = ^^ ⋅ cos ^^^^sin

[0119] Assuming the walk ^^^^has the Markovproperty, since: ^^(^^^^|^^^^−1) = ^^(^^^^|^^^^−1, ^^^^−2, ... , ^^1)

[0120] As long as we pick a stopping time ^^ for ^^ that is a multiple of ^^, if we define a new process ^^(^^) ^^as: ^^(^^)^^ = ^^^^+^^ − ^^^^

[0121] we will have the strong Markov property for ^^(^^) ^^, and so it will be independent of ^^^^. with the steps ^^^^distributed identically to thethe ^^^^random walk.

[0122] We also have the density of the ^^^^for a^^ is a normal distribution by CLT. Moreover, since the distribution of the ^^^^has mean zero, we will show the sample paths of the ^^^^are self-similar by the strong Markov property, which allows us to use the reflection principle to derive a density for the hitting times. Much of the following work mirrors the proof of the reflection principle in Brownian Motion, Martingales and Stochastic Calculus by Jean-Francois Le Gall. (Gall, J. L. (2016), Brownian Motion, Martingales, and stochastic calculus, In Graduate texts in mathematics, https: / / doi.org / 10.1007 / 978-3-319-31089-3).

[0123] Using the strong Markov property of ^^, we can prove the following: For every^^ ∈ ℕ, set ^^^^^^ = max^^^^≤^^^^^^^^^^. Then if ^^ ≥ 0 and ^^ ∈ (−∞, ^^] we have^^(^^^^^^ ≥ ^^, ^^^^^^ ≤ ^^) = ^^(^^^^^^ ≥ 2^^ − ^^))20 Docket No. JTS0005

[0124] And ^^^^^^has the same distribution as |^^^^^^|. Then, by Chung-Fuchs theorem,we have that ^^^^ < ∞, where ^^^^ is the first hitting time of ^^^^ for distance ^^. Then we have:^^(^^ ≥ ^^, ^^ ≤ ^^) = ^^(^^ ) (^^^^)^^^^ ^^^^ ^^ ≤ ^^^^, ^^^^^^ = ^^ (^^^^ ≤ ^^^^, ^^^^^^−^^^^ ≤ ^^ − ^^^^^^)

[0125] ^^ (−^^ (^^ ) ≤ ^^ − ^^ ) = ^ ( ^^)^^^^−^^^^ ^^^^ ^ (^^^^^^−^^^^ ≤ ^^ − ^^^^^^)

[0126] Then, ^^}, we can^^ (^^^^ ≤ ^^, −^^ (^^^^^^^^−^^^^ ≤ ^^ − ^^^^^^) = ^^(^^^^ ≤ ^^, ^^^^^^ ≥ 2^^^^^^ − ^^) = ^^(^^^^^^ ≥ 2^^^^^^ − ^^)^^^^. Note, in the Brownian case, since uncountably many random Gaussian variablesgenerate the process, ^^^^^^ = ^^ is justified by the continuity of sample paths. Here, sincerandom vector defining the process ^^^^ given by (^^1, ^^2, ... ^^^^) is countable, we cannot saythat ^^ = ^^, but that ^^^^^^ ≥ ^^. However,variation of the ^^ is bounded,can say that:^^(^^^^^^ − ^^ > ^^) < ^^,

[0128] for some small ^^, ^^.^^(^^^^^^ ≥ 2^^^^^^ − ^^) ≈ ^^(^^^^^^ ≥ 2^^ − ^^)

[0129] Then we^^(^^^^^^ ≥ ^^) = ^^(^^^^^^ ≥ ^^, ^^^^^^ ≥ ^^) + ^^(^^^^^^ ≥ ^^, ^^^^^^ ≤ ^^)^^ ^^ ^^ ^^ ^^

[0130] whereand {^^^^^^≥^^} ⊂ {^^ ≥ ^^}. We can now derive t ^ ∑^^ ^^^^ he density of the ^^^ . We first characterize ^^=1 ^^^^to define the variance of the ^^ . We have that for each ^^^^, each step is takenthe update equations above, where we takesteps ^^ corresponding to a small distance along the path of ^^. These small steps in distance correspond to a small step intime, according to the velocity of light in the medium. Then, given ^^ ≪ 1, we have that 1 / ^^is the number of steps taken at time ^^ = 1 (note in our calculation of the ^^^^, our choice of^^ impacts ^^1 ^^.). Then we have ^^ =^^^^is the number of ^^^^ steps we take by time ^^ = 1. Then21 Docket No. JTS0005we have that ∑^^ ^^=1 ^^^^will be approximately normally distributed by CLT and will have some variance ^^2.

[0131] Then,the property of normal distributions that samples ^^^^ from ^^(0, ^^2)satisfy ∑^^ ^^=1 ^^^^ ∼ ^^(0, ^^^^2), and √^^^^^^ ∼ ^^(0, ^^^^2), we have that ^^^^has the same ^^

[0132] the equation above, we have:^^ ≤ ^^^^) = ^^(^^^^^^ ≥ ^^) = ^^(|^^^^^^| ≥ ^^) = ^^(^^ 2 2^^^^ ≥ ^^ ) =^^^^ ^^2^^^^

[0133] Then, ^^2^^^^^^ (^^^2^≤ ^^ ).We consider ^^(^^) = ^^2, and so ^^−1(^^)√^^2. Since ^^ is decreasing on its^^^2^^^^^−1 is likewise{^^(^^) ≤^^^^ ^^} = {^^ ≥ ^^−1 (^^^^ ^^ )} and so we have^^ (^^(^^) ≤^^^^) = ^^ (^^ ≥ ^^−1 (^^^^ )).^^ ^^ via^^ ,20 ^^1^^2^^^^

[0135] Butexpression by^^2multiplying by 2, which accounts for the concentration of probability mass due to^^2mapping all values of ^^ into [0, ∞).

[0136] Thus, we have that the density of hitting times ^^^^is given by: ^^ ^^ = ^^ =^^ ^^−3 / 2 ^^2 / (2^^^^2)

[0137] Where we^^

[0138] FIG. 7 provides empirical CDF of samples generated from Monte-Carlo 22 Docket No. JTS0005simulation of the photon update equations compared to the modified Levy distribution, using ^^ learned from the MC data with a PyMC Bayesian sampler.

[0139] To test the analysis above, we generate samples from ^^^^(^^|^^, ^^),Bernoulli(^^^^), and Uniform(0,2^^), and simulate 20,000 photons for 5 × 108steps, simulated using a photon scattering algorithm. The algorithm was implemented in Jax, using the vmap parallelization functionality and just-in-time compilation features.20,000photon packets were simulated for 5 × 108 steps in 20 minutes on an RTX-3090.

[0140] Here, each photon consists of an ^^, ^^, ^^ coordinate, a current heading, a photonstrength value, and a hitting time value. We start each photon with a strength of 100, and every time a photon is absorbed, we subtract 1 from the strength. Photon positions are updated according to the update equations referenced above.

[0141] We then calculate the cumulative probability distribution of the hitting times at^^ = 0 for the Monte-Carlo particles and compare it to the CDF of the modified Levydistribution with parameters ^^ = 1, and ^^ = 1.29 learned by Bayesian estimation from theMonte-Carlo samples. The two are plotted in FIG.7.

[0142] We note here the uptick in samples at the very end. This is due to the simulationtime being capped at 5 × 108 steps. If we were to simulate for longer, it would bereasonable to expect the tails to match up correctly, given a large enough number of particles. Here, we note the maximum error between this final value of the MC samples and the Levy distribution is 0.009 or about 17%. This error is should be acceptable for most purposes, as larger times mean more attenuation, so error from tail photons is negligible.

[0143] These results support the assumption that the modified Levy distribution describes the hitting times ^^^^, which correspond precisely to the ^^^^in the equation above. Thus, we have a simple summation to represent the incident light at any particular wavelength, and since we can describe ^^^^and ^^^^with probability density, we can further approximate the summation to give a closed-form model for incident light.

[0144] It may be desirable to further provide an integrated recursive path approximation. If we wish to approximate the sum [IRPA_SUM], using two terms, and account for the time-varying nature of the signal, due to the pulse of the blood. We wish 23 Docket No. JTS0005to split up the sum based on path length, and after more analysis, we can show: ^^^^on t,since we now consider the intensity across time. ^^^^has no dependence on time since ^^^^determines how long a photon path spends in the body, not what it encounters. ^^^^accounts for this change in path absorption across time.

[0146] We have that ^^1, ^^2 sum to 1 and are the integrals of the modified Levydistribution from 0 to some time ^^^^corresponding to paths that have a near zero probability of entering the artery, and the integral of the modified Levy from ^^^^to ∞, respectively. Thus, we have: ^^^^^^1 = ∫ ℒ (^^|^^, ^^)^^^^

[0147] Then ^^1, ^^2 are thedetected path over the twosections of the Levy distributions, respectively. Thus, we have: ^^−^^^^1^^^^^^11^^^^=^^−^^^^^^^^^^1 ℒ ^, ^^)^^^^

[0148] where ^^ 1 , ^^^^^^2path bin,Note we scale by 1 / ^^1, 1 / ^^2 respectively so each probability density measure over thesubset of path lengths still integrates to 1.

[0149] Finally, ^^(^^)is some function that determines what percentage of time the ^^2path spends in tissue vs blood, dependent on the pulse of the blood, which is dependent on time. 24 Docket No. JTS0005

[0150] Note ^^^^ and ^^(^^) are codependent since ^^^^ determines ^^1, ^^2, ^^1, ^^2. To expandthe analysis on ^^(^^), we can decompose it as ^^(^^) = ^^ ⋅ ^^(^^) + ^^, where ^^(^^) is strictly afunction of the pulse, which models the volumetric change ofof the body we monitor. We have that ^^ is an offset modeling the percentage of time paths spend in the tissue at both systolic and diastolic blood volumes, and ^^ is a scalar that parameterizes the change in percentage of time paths spend in the blood between systolic and diastolic blood volumes.

[0151] Here, we note that ^^ is a proxy for the distance of the artery from the skin. Then, given some parameter ^̂^ (which might have to be learned) we have for the distanceparameter ^^ in the Levy distribution: ^^ = ^̂^ ⋅ ^^. We can also define ^^ as an offsetparameter so that paths start some distance ^^ into the skin. This is done so the Levy distribution does not collapse to a Dirac distribution. For most purposes, it is sufficient to set ^^ to 1 (assuming a length unit of millimeters).

[0152] Then, given ^^, we may determine an upper bound on the probability that a path of length ^^ enters the artery. We define this function ℒ: ^^1^^2ℒ(^^|^^, ^^, ^^) = sup {∫ ℒ (^^|^^, ^^) ∫ ℒ (^^|^^ + ^^, ^^)}^^1,^^2|

[0153] Note here this expressionbut we wish toconsider not ^^^^ only at ^^ = ^^^^, but inf^^≤^^^^{^^^^}, as we could have ^^^^ enter the artery atsome time ^^ < ^^^^, but then rise back above the level of the artery. By self-similarity wecan consider the suprema of two shorter paths, which leads to our definition of ℒ. To justify this definition, we can consider the diagram of FIG.8.

[0154] In FIG.8, the ^^1corresponds to the ^^1in the upper bound of integration on the leftmost integral in the definition of ℒ, andlikewise on the rightmost integral in the definition of ℒ.

[0155] Then ℒ(^^) computes the maximal probability that a path both enters the artery and returns to the skin, given that the total path length must be less than ^^. We consider two sub-paths, where the first path represents the distance from the original path’s starting point, and the second is the distance from the artery to the skin. This reduces to evaluating the density of hitting times of paths of distance ^^ corresponding to the ^^1path 25 Docket No. JTS0005and the density of hitting times of paths of distance ^^ + ^^, corresponding to the ^^2 path.Then, the product of these two integrals gives the probability that a path with length less than or equal to ^^1hits distance ^^ and that a path of length less than or equal to ^^2hitsdistance ^^ + ^^.

[0156] Then, because the paths exhibit self-similarity and the strong Markov property, we have independence between time steps, so we can factor the probability that a path with length less than or equal to s both hits the artery and the skin as the product of the probability that two paths hit the artery, and the skin, which leads to the expression involving the product integrals over the Levy density. This gives ℒ(^^) as the upper bound on this probability for a given path length ^^.

[0157] Then once ^^ is determined, we have that given ^^, we may compute ^^^^as: ^^^^ = sup {^^|ℒ(^^|^^, ^^, ^^) < ^^}^^

[0158] Where ^^ is some the probability that a path oflength ^^^^has at most probability ^^ of entering the artery and returning to the skin.

[0159] Then ^^^^ is determined given ^^, and once ^^^^ is known, then ^^1, ^^2, ^^1, ^^2 aredetermined, which necessarily determines ^^ in the expression ^^(^^) = ^^ ⋅ ^^(^^) + ^^, given^^(^^).

[0160] Preferred embodiments and implementations may utilize machine learning models that exploit this new understanding to extract parameters of interest, such as ^^^^.This is because ^^^^ is a weighted average of ^^^^, ^^ℎ..., which are the extinction coefficientsof red blood cells and other blood parameters, and the weights on those coefficients correspond to the absolute concentration of those parameters of interest. In some cases, the common term ^^0^^1^^−^^^^1^^^^may be isolated, which might be done using the loss function: ^^log(^^^^1^^ − ^^) / ^^^^log(^^^^2^ − ^^) / ^^argmin^^^^ ^^^^^^^,^^ var (^^) + var () log(^^ − ^^^^ ^^^^^^2^^ ) / ^^^^^^log(^^^^1^^ − ^^) / ^^

[0161] where ^^^^1^^ , ^^^^2^^ correspond to the detected luminance at time t for wavelengths1 and 2, respectively. Here var is the sample variance function, which should be minimized when x and y correspond to ^^−^^^^0 1^^^^11^^1^^, ^^02^^1^^−^^^^1^^^^2, respectively. This is 26 Docket No. JTS0005because the fraction: ^^log(^^ − ^^ ^^ ^−^^^^1^^^^1 −^^^^1^^^^1^^^^^^1^^ 01 1 ^) / ^^01^^1^^ ^^

[0162] evaluates to a short period of time, andso should have zero variance over a short time (one or two heartbeats). One might replace the sample variance by the information-theoretic differential entropy function, as a constant value will have zero entropy (Dirac distribution).

[0163] Experimental support for the integrated model disclosed herein has been established. Initial results on real data using the loss function defined above are described below and are shown in FIGS.9A and 9B. We gathered real world spectrometer data in order to test the effectiveness of the integrated model developed above.

[0164] The experimental data ^^^^1,^^ , ^^^^1,^^ are taken from an Hamamatsu InGaAs (IndiumGalium Arsenide) infrareddeveloped by Tula Health, Inc. in conjunction with the Brigham Young University Physics Department. To produce the data used below, the spectrometer was placed against the radial artery in the wrist of a test subject, and data were recorded to facilitate further analysis.

[0165] Given this data, we choose ^^^^1at wavelength 1300 nm, and ^^^^2at wavelength1100 nm. In FIG. 9A, we see the loss plotted for various values of ^^, ^^ corresponding tovarious DC offset values.

[0166] In FIG.9B, we see a similar loss surface generated for synthetic data generated using the equations and techniques described herein. Here, the absolute magnitudes of the loss surface are irrelevant, as the synthetic data have less noise than the experimental data, which increases the variance for all loss values, so we only consider the relative magnitudes of the real vs. synthetic data.

[0167] In FIGS.9A and 9B, it can be seen that the loss surfaces of the experimental and synthetic data are qualitatively similar, which is an indication that the integrated model captures much of the structure of experimental data. This suggests that the integrated model is representative of much of the structure seen in time dependent photon scattering measurements. 27 Docket No. JTS0005

[0168] Lemma 2:

[0169] The sum of two uniform random variables supported on the unit circle will be uniformly distributed.

[0170] Proof:

[0171] Because circular convolution is equivalent to identifying the endpoints of the convolution interval, we may simply compute the circular convolution of two uniform random variables.

[0172] Thus, let ^^ and ^^ be independent random variables, each uniformly distributed on the interval [0,1). The circular convolution of these two distributions is defined as thedistribution of the random variable ^^ = ^^ + ^^ mod 1.

[0173] The probability density function (pdf) of a uniform distribution on [0,1) is given by: ^^^^(^^) = ^^^^(^^) = {1 if 0 ≤ ^^ < 1,0 otherwise.

[0174] The pdf of the sum ^^ = ^^ + ^^ (before taking the modulo operation) is given bythe convolution of ^^^^and ^^^^: ∞ ^^^^(^^) = ∫ ^^^^ (^^)^^^^(^^ − ^^) ^^^^.−∞

[0175] Since ^^^^(^^) = ^^^^(^^) = 1 for ^^ ∈ [0,1) and 0 otherwise, the convolution integralsimplifies to: min(^^,1)^^^^(^^) = 1  ^^^^ = min 1) − max ^^ − 1).

[0176] Now,the pdf of ^^ mod 1, which we will denote as ^^^^%1(^^). The pdf ^^^^%1(^^)is periodic with period 1, so we only need to consider ^^ in the interval [0,1):

[0177] For 0 ≤ ^^ < 1, ^^^^%1(^^) = ^^^^(^^) + ^^^^(^^ + 1). In this interval, ^^^^(^^) = ^^ and^^^^(^^ + 1) = 0, so ^^^^%1(^^) = ^^ + 0 = ^^.

[0178] For 1 ≤ ^^ < 2, ^^^^%1(^^ − 1) = ^^^^(^^ − 1) + ^^^^(^^). In this interval, ^^^^(^^ − 1) = 0and ^^^^(^^) = 1 − (^^ − 1) = 2 − ^^, so ^^^^%1(^^ − 1) = 0 + (2 − ^^) = 2 − ^^.

[0179] Combining these two cases, we find that the pdf of ^^ mod 1 is constant andequal to 1 for ^^ ∈ [0,1), which is the pdf of a uniform distribution on [0,1). Therefore, the28 Docket No. JTS0005circular convolution of two uniform distributions is also uniform.

[0180] Input: Boundary conditions, absorption and scattering probabilities, path length^^, time step ^^^^, scattering parameters ^^, ^^ Output: Updated particle state, final particle.

[0181] Sample action from categorical distribution Determine new direction using particle position = particle position = subtract 1 from photon strength particle position =

[0182] Generate sampling ^^1, ^^2 from Uniform(0,1),Uniform(0,2^^) Find ^^ bycomputing ^^^^−1(^^1|^^, ^^) Update new heading by computing ^^^^^^^^, ^^^^^^^^ from^^cos(^^2), ^^sin(^^) new direction

[0183] compute direction ^^, ^^ as current particle heading + new heading compute newposition as ^^^^^^^^ = ^^^^^^^^ + ^^sin(^^)cos(^^), ^^^^^^^^ = ^^^^^^^^ + ^^sin(^^)sin(^^), ^^^^^^^^ = ^^^^^^^^ + ^^cos(^^)^^^^^^^^ , ^^^^^^^^ , ^^^^^^^^

[0184] FIGS.10A and 10B illustrate schematically a system for using electromagnetic radiation from an electromagnetic emitter 1010 to process signals and track various physiological conditions. As shown in FIG. 10B, radiation is emitted from emitter 1010 and is directed through the skin 1030 of a user, through various tissue layers 1035, and into a blood vessel 1040. The radiation is reflected at various points and received by a sensor 1020, which may comprise, for example, a spectrometer. Thus, the signal data referenced throughout this disclosure may be received by a physiological sensor, such as a non-invasive blood monitor. Examples of such monitors can be found in U.S. Patent No.11,903,686 titled “Systems, Apparatuses, and Methods for Optimizing a Physiological Measurement Taken From a Subject,” the entire contents of which are hereby incorporated herein by reference in its entirety.

[0185] As shown in FIG. 10B, the sensor 1020 may be configured to measure the magnitude / intensity of the received radiation over time. The resulting signal may therefore consist of a pulsatile component 1002, which corresponds to radiation 1002 being reflected from blood vessel 1040, as shown in FIG.10B. Because the blood vessel 1040 contracts and expands at a rate corresponding with the heartbeat of the user, this portion of the signal may contain a periodic and / or pulsatile signal, as indicated by the sine wave within section 1002 of the signal.

[0186] The signal may also include a non-pulsatile component 1004 that results from 29 Docket No. JTS0005radiation reflected from within the blood vessel 1040. In systems configured to predict, estimate, and / or detect blood analytes, non-pulsatile blood component 1004, otherwise referred to herein as the “baseline” component, is likely to contain useful information that may be used to provide a better prediction or estimation of the blood analyte under consideration.

[0187] The signal may further comprise a DC offset or tissue component 1006. DC offset / tissue component 1006 is the result of radiation being reflected from non-blood tissue 1035 of the user. Again, if the desired result is to assess blood analyte conditions, this portion of the signal is not likely to contain useful information. Moreover, having this portion of the signal present may detract from the useability of components 1002 and 1004. This is because prior techniques for processing data with all three components 1002 / 1004 / 1006 fail to differentiate between components 1004 and 1006. For example, it has been typical to simply assume that components 1004 and 1006 were both resulting from tissue reflections, and therefore did not contain useful information. However, the techniques provided herein allow for not only differentiation of components 1004 and 1006, but also for extraction of one or both of these components. For example, it may be desirable to extract DC offset / tissue component 1006 entirely so that it does not detract from the usefulness of the information in the other two components 1002 and 1004.

[0188] FIG.11 is a flowchart illustrating an example of a method 1100 for processing signal data to estimate, calculate, and / or predict a condition, such as a blood analyte condition, according to some implementations. At step 1102, signal data may be received. In some cases, the signal data may be multi-channel signal data. For example, each channel of such a multi-channel signal may correspond with a distinct wavelength, wavelength range, and / or sensor. This may provide for enhanced detection capabilities by taking advantage of different absorption characteristics of glucose or other analytes at various wavelengths, typically within the near-infrared (NIR) spectrum. As previously mentioned, the signal data may be received from a blood monitoring device and / or system, which may comprise a wearable device, such as a wrist monitor, that may be configured with the aforementioned emitter(s) and sensor(s).

[0189] Although use of a single channel of signal data may be feasible for some uses, a multi-channel signal may allow for capturing the signal Si(t) at different settings. 30 Docket No. JTS0005

[0190] According to the preferred embodiments, these signals may have the functional form: Si(t) = Ci + Di · exp(−vi · p(t)), where Ci is the offset value for each channel I, Di is an unknown scaling factor, vi is a vector of weights dependent on the channel i. and p(t) is a pulse function representing the physiological signal.

[0191] By collecting multi-channel data, multiple observations of the same underlying physiological process may be accessed and / or may be modulated differently according to the weights vi. This setup allows leveraging the relationships between channels to extract the common offset C and improve signal processing.

[0192] At step 1104, the incoming signals may be processed to remove the mean offset. In some implementations and embodiments, the offset C may be removed from each signal by subtracting a proposed offset value c. Since C is unknown, in some cases, one could start with an initial guess c and adjust it through optimization. After subtracting c, the natural logarithm of the result may be taken. According to the preferred implementations of the model, this yields: log(Si(t) − c) = log(D) − vi · p(t).

[0193] Since D is unknown and acts as an additive constant in the logarithmic domain, in some cases it can be ignored, instead focusing on frequency components greater than zero.

[0194] Subtracting c aims to estimate and remove the offset C. Taking the logarithm linearizes the exponential function in the model, transforming the multiplicative relationship into an additive one. Ignoring log(D) (the DC component) in subsequent frequency-based analysis allows us to compare the dynamic parts of the signals across channels.

[0195] The signals may then be transformed to the frequency domain at step 1106. In some implementations, the logarithmically transformed signals may be converted into the frequency domain using the Fast Fourier Transform (FFT). This transformation enables the analysis of frequency components associated with the pulsatile function p(t).

[0196] Although not required for all contemplated embodiments and implementations, transforming to the frequency domain isolates the frequencies related to p(t). Removing the DC component (log(D)) ensures focus on the variations caused by −vi * p(t), which differ across channels due to different vi values.

[0197] One or more features, parameters, and / or components may then be extracted 31 Docket No. JTS0005from the data at step 1108. In some cases, this may be done via machine learning. In some implementations, step 1108 may comprise excluding the zero frequency (DC) component and / or the highest frequency component, which may allow for concentrating on the frequencies that carry the pulsatile information. This step may also isolate the relevant features influenced by p(t) and vi.

[0198] By removing the zero-frequency component, one can effectively eliminate log(D), a constant offset. Focusing on non-zero frequencies allows for analyzing the dynamic behavior of the signals due to p(t).

[0199] The signal(s) may then be adjusted, aligned, and / or normalized at step 1110. In some cases, step 1110 may comprise reducing scale discrepancies in the signal(s). In some cases, the signals may be normalized and aligned to account for sign differences in vi. Since vi may have different signs across channels, aligning ensures that the signals are comparable. Normalization may remove differences in scale, allowing direct comparison of signal shapes.

[0200] Alignment may be useful for correcting discrepancies due to different signs of vi. Normalizing the signals removes scaling differences, ensuring that any remaining variance is due to incorrect offset c rather than amplitude differences.

[0201] Step 1112 may comprise computing discrepancy measure and / or variance across the signal channels. After alignment and normalization, the variance across channels may be computed, preferably at each frequency. If c = C, the transformed signals should be identical across channels, resulting in zero variance. Nonzero variance indicates that the offset c is incorrect.

[0202] Because the variance serves as a discrepancy measure or loss function, minimizing this variance helps find the value of c that best estimates the true offset C, thereby leveraging the model’s relationships between channels.

[0203] The offset adjustments may then be improved and / or optimized to reduce and / or minimize discrepancy and / or variance at step 1114. In some implementations and related embodiments, the computed variance (loss function) may be used in an optimization routine to adjust the proposed offset c. By iteratively updating c to minimize the loss, convergence towards the true offset C may be achieved.

[0204] Optimization techniques (e.g., gradient descent) may use the loss to update c. 32 Docket No. JTS0005Since the loss is differentiable with respect to c, gradients can be computed and the c that minimizes the variance, or at least reduces the variance, can be identified, thereby allowing for more accurately estimating C.

[0205] Method 1100 may conclude when the optimization converges and the variance reaches a minimum value, or a predetermined minimum threshold value.

[0206] FIG.12 is a schematic diagram of a system 1200 for processing data, such as data from a wearable device 1212 used for monitoring health conditions, such as blood analyte conditions, according to some implementations. As shown in this figure, wearable device 1212, which may comprise a wrist monitor in some embodiments, comprises one or more electromagnetic emitters 1214 and one or more sensors 1216 configured to receive reflected radiation from the emitter(s) 1214. In some embodiments, the wearable device 1212 may be positioned on the wrist so that a physiological sensor 1216 is positioned adjacent to a particular, targeted physiological feature, such as an ulnar artery or radial artery, for example.

[0207] In some embodiments, an electromagnetic emitter 1214 may be positioned at a suitable location relative to the sensor1216 so that electromagnetic radiation of a desired wavelength may be emitted from emitter 1214, reflected from desired regions within the body, such as along an adjacent to a particular artery in the case of blood monitoring, and then received by sensor 1216 for processing according to, for example, the method of FIG.11 or any of the other methods disclosed herein.

[0208] As those of ordinary skill in the art will appreciate, sensor 1216 may be configured to convert reflected / incoming radiation into electronic signals, which may be transmitted to and received by processing system 1201. For example, in some embodiments, sensor 1216 may be configured to measure the intensity of the received radiation over time and transmit this signal to system 1201. System 1201 may be incorporated onto the wearable device 1212 or may be incorporated into another device, such as a mobile application running on a smartphone, a computer, a remote server, or the like.

[0209] System 1201 may comprise one or more processors 1202. Processor(s) 1202 may be configured to receive and / or process data from wearable device 1212 and / or sensor 1216. In addition, processor(s) 1202 may be configured to perform various 33 Docket No. JTS0005functions, some of which may be part of a particular software module. As used herein, a software module or component may include any type of computer instruction or computer executable code located within a memory device and / or m-readable storage medium. A software module may, for instance, comprise one or more physical or logical blocks of computer instructions, which may be organized as a routine, program, object, component, data structure, etc., that perform one or more tasks or implements particular abstract data types.

[0210] System 1201 may comprise a Fourier Transform module 1204. Fourier Transform module 1204 may be configured to take raw and / or preprocessed data, which preferably includes a periodic feature, and project the signal data onto a set of basis functions. In some embodiments, Fourier Transform module 1204 may be configured to extract one or more features of a periodic aspect of the signal, such as features associated with a shape of a pulsatile signal contained within the data. In some embodiments, certain frequency amplitudes computed by the Fourier Transform Module 1204 might be indicative of certain analyte concentrations, but also may be indicative of certain cardiovascular conditions. It may therefore be possible to diagnose, for example, atrial fibrillation or similar conditions by analyzing the shape of the waveform with Fourier analysis.

[0211] System 1201 may further comprise a Tissue-Dependent Signal Extraction module 1206. In some embodiments, Tissue-Dependent Signal Extraction module 1206 may be configured to utilize one or more of the steps described above in connection with method 1100 to extract a DC offset portion of a signal, which in some cases may be the portion of the signal corresponding with reflections from non-blood tissue, or may otherwise comprise portions of the signal that may not be needed and / or may be detrimental to providing accurate assessments of blood analyte content. In some embodiments, Tissue-Dependent Signal Extraction module 1206 may therefore be configured to use one or more of the DC extraction methods described herein.

[0212] System 1201 may further comprise a Non-Pulsatile Blood Signal Extraction module 1208. In some embodiments, module 1208 may be configured to extract a portion of the signal corresponding to reflections from the blood itself rather than the blood vessel wall. In some cases, module 1208 may comprise simply compiling the data remaining 34 Docket No. JTS0005following extraction of the DC offset portion of the signal, along with the pulsatile signal (described below).

[0213] Some embodiments may further comprise a Pulsatile Signal Extraction module 1210. In such embodiments, module 1210 may be configured to extract the portion of the signal that comprises a periodic and / or pulsatile signal received by sensor 1216. In embodiments in which system 1200 is used to provide blood analyte condition data, such as glucose monitoring data, module 1210 may therefore be configured to extract the portion of the signal associated with a heartbeat signal or other physiological, periodic signal. As those of ordinary skill in the art will appreciate, in some embodiments, module 1210 may utilize Fourier transforms and / or analyses, such as the Fast Fourier Transform (FFT). In some embodiments and related implementations, once the DC offset has been extracted, as described throughout this disclosure, the remaining signal, which comprises a pulsatile signal and a non-pulsatile portion, which may comprise a non-pulsatile blood portion in some cases, may be processed to separate them by simply extracting the offset from the remaining signal and / or extracting the pulsatile portion, either of which may be accomplished in a relatively straightforward manner.

[0214] Another method 1300 according to some implementations is illustrated in the flowchart of FIG. 13. In this method, signal data is received at step 1302. In some implementations, the signal data may comprise multi-channel signal data. In some implementations, the signal data may be received from a physiological sensor and / or monitor, such as a blood glucose monitor or a pulse oximeter, for example.

[0215] The signal data may then be processed at step 1304 to extract one or more components. For example, in some implementations in which blood analytes, such as oxygen or glucose, are being monitored, the signal may comprise useful data and non- useful data. If the device that is providing the data signal is configured to use electromagnetic radiation directed towards a blood vessel, for example, it may be useful to extract a component, or at least as much of a component as possible, that corresponds to reflected radiation from tissue rather than blood and / or blood vessels. This component may otherwise be referred to herein as a DC offset component. Techniques for extracting this component are taught throughout this disclosure.

[0216] In addition, for some purposes, it may be desirable to extract a non-pulsatile 35 Docket No. JTS0005blood component from the signal. As previously mentioned, this component may correspond with reflected radiation from the blood itself within a blood vessel. Because techniques for extracting the pulsatile signal and the DC offset / tissue component are taught herein, the remaining portion of the signal may be considered the non-pulsatile blood component and may therefore be extracted by extracting the other two components. It is thought that this component may be used alone for certain purposes.

[0217] A condition, such as a blood analyte condition, may then be calculated, predicted, estimated, and / or assessed at step 1306. In some implementations, this step may be performed by use of machine learning algorithms. Additionally, or alternatively, this step may be performed by calibrating the processed and / or extracted data, such as by calibrating results with more accurate and / or invasive techniques for assessing blood analyte conditions. In some cases, calibration may be achieved by finding blood analytes using known lab procedures, such as metabolic panels and the like, which may then be compared with the results of the modeling and / or processed data.

[0218] In certain embodiments, a particular software module may comprise disparate instructions stored in various locations of a memory device, which together implement the described functionality of the module. Indeed, a module may comprise a single instruction or many instructions, and may be distributed over several different code segments, among different programs, and across several memory devices. Some embodiments may be practiced in a distributed computing environment where tasks are performed by a remote processing device linked through a communications network. In a distributed computing environment, software modules may be located in local and / or remote memory storage devices.

[0219] In addition, data being tied or rendered together in a database record may be resident in the same memory device, or across several memory devices, and may be linked together in fields of a record in a database across a network. Furthermore, embodiments and implementations of the inventions disclosed herein may include various steps, which may be embodied in machine-executable instructions to be executed by a general-purpose or special-purpose computer (or another electronic device). Alternatively, the steps may be performed by hardware components that include specific logic for performing the steps, or by a combination of hardware, software, and / or firmware. 36 Docket No. JTS0005

[0220] Embodiments and / or implementations may also be provided as a computer program product including a machine-readable storage medium having stored instructions thereon that may be used to program a computer (or other electronic device) to perform processes described herein. The machine-readable storage medium may include, but is not limited to: hard drives, floppy diskettes, optical disks, CD-ROMs, DVD- ROMs, ROMs, RAMs, EPROMs, EEPROMs, magnetic or optical cards, solid-state memory devices, or other types of medium / machine-readable medium suitable for storing electronic instructions. Memory and / or datastores may also be provided, which may comprise, in some cases, non-transitory machine-readable storage media containing executable program instructions configured for execution by a processor, controller / control unit, or the like.

[0221] It will be understood by those having ordinary skill in the art that changes may be made to the details of the above-described embodiments without departing from the underlying principles presented herein. Any suitable combination of various embodiments, or the features thereof, is contemplated.

[0222] Any methods disclosed herein comprise one or more steps or actions for performing the described method. The method steps and / or actions may be interchanged with one another. In other words, unless a specific order of steps or actions is required for proper operation of the embodiment, the order and / or use of specific steps and / or actions may be modified.

[0223] Throughout this specification, any reference to “one embodiment,” “an embodiment,” or “the embodiment” means that a particular feature, structure, or characteristic described in connection with that embodiment is included in at least one embodiment. Thus, the quoted phrases, or variations thereof, as recited throughout this specification are not necessarily all referring to the same embodiment.

[0224] Similarly, it should be appreciated that in the above description of embodiments, various features are sometimes grouped together in a single embodiment, figure, or description thereof for the purpose of streamlining the disclosure. This method of disclosure, however, is not to be interpreted as reflecting an intention that any claim require more features than those expressly recited in that claim. Rather, inventive aspects lie in a combination of fewer than all features of any single foregoing disclosed 37 Docket No. JTS0005embodiment. It will be apparent to those having skill in the art that changes may be made to the details of the above-described embodiments without departing from the underlying principles set forth herein.

[0225] Likewise, benefits, other advantages, and solutions to problems have been described above with regard to various embodiments. However, benefits, advantages, solutions to problems, and any element(s) that may cause any benefit, advantage, or solution to occur or become more pronounced are not to be construed as a critical, a required, or an essential feature or element. The scope of the present invention should, therefore, be determined only by the following claims.

Claims

embodiment. It will be apparent to those having skill in the art that changes may be made to the details of the above-described embodiments without departing from the underlying principles set forth herein. [00225] Likewise, benefits, other advantages, and solutions to problems have been described above with regard to various embodiments. However, benefits, advantages, solutions to problems, and any element(s) that may cause any benefit, advantage, or solution to occur or become more pronounced are not to be construed as a critical, a required, or an essential feature or element. The scope of the present invention should, therefore, be determined only by the following claims. 38 Docket No. JTS0005CLAIMS 1. A method for processing signals from a non-invasive blood analyte tracking system, the method comprising the steps of: receiving signal data from a non-invasive blood monitor, wherein the signal data comprises: a frequency component corresponding to a heart rate pulse of a user of the non-invasive blood monitor; a non-pulsatile blood component corresponding to light reflected from within a blood vessel; and a tissue-dependent DC offset component; processing the signal data to extract the tissue-dependent DC offset component; and evaluating the processed signal data to provide an estimate of a blood analyte condition using only the frequency component and the non- pulsative blood component from the signal data.

2. The method of claim 1, wherein the blood analyte condition comprises a concentration of the blood analyte.

3. The method of claim 1, wherein the blood analyte condition comprises a trend over time associated with the blood analyte.

4. The method of claim 1, wherein the blood analyte comprises glucose.

5. The method of claim 1, wherein the step of processing the signal data comprises use of a probability density model.

6. The method of claim 5, wherein the probability density model comprises a Levy Distribution. 39 Docket No. JTS00057. The method of claim 1, wherein the blood analyte comprises oxygen.

8. A method for extracting a DC offset component from signal data using a non-invasive blood monitor, the method comprising the steps of: receiving signal data from a non-invasive blood monitor; processing the signal data to remove a mean offset component from the signal data; transforming the processed signal data to a frequency domain; extracting a DC offset component from the transformed signal data; and using unextracted data to predict a blood analyte condition.

9. The method of claim 8, wherein the blood analyte condition comprises a concentration of the blood analyte.

10. The method of claim 8, wherein the blood analyte condition comprises a trend over time associated with the blood analyte.

11. The method of claim 8, wherein the blood analyte comprises glucose.

12. The method of claim 8, further comprising adjusting signal data to reduce scale discrepancies.

13. The method of claim 8, further comprising processing the signal data to compute discrepancy measure.

14. The method of claim 8, further comprising adjusting the offset to reduce discrepancy. 40 Docket No. JTS0005

Images