% SIMIODETemplate-TeacherORStudent.tex and simiode.cls for LaTeX2e
% Version 29 July 2015.
%
% This class (simiode.cls) is a modification of that used by
% the journal PRIMUS - Problems, Resources, and Issues in
% Mathematics Undergraduate Studies.
%
% See http://www.tandfonline.com/loi/upri .
%
% The class was developed by  Marc Frantz with suggestions of the
% PRIMUS Editorial Board, especially John Thoo.
%
% This template is to be used with the  class file "simiode.cls".
%
% The latest version of this class and template files
% can be found at
%
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
%
\documentclass{simiode}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Basically the Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review (but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.
% Rename the pdf after you compile each final version.
% This option was created  by Eric Sullivan, Carroll College, Helena MT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will
% look in final form, so that you can better plan your
% presentation.  To see your manuscript formatted very
% similarly to the on-line published format of  SIMIODE,
% comment the following line out by inserting a percent sign % (%) at its beginning.

\renewcommand{\baselinestretch}{1.3}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding macros %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Any macros defined for your paper should be contained in
% the top matter. Likewise, any environment definitions such
% as \newtheorem or \newenvironment should also go in the
% top matter.
%
% Do not \input your macros as separate files within your
% final version because more files make it harder for the
% editors and users to keep track of and/or modify and
% append your paper.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\newpage


\title{\ifnum\TeacherEdition=1 TEACHER VERSION  \else STUDENT VERSION \fi\\MODELING INTRAOCULAR GAS BUBBLES IN\\RETINAL SURGERY PATIENTS}



\markboth{Eye Model}{Eye Model}





\author{Brian Winkel,
Director SIMIODE\\
Chardon OH 44024 USA\\
BrianWinkel@simiode.org
}
\makeStitlePDFLaTex


\begin{abstract}
We offer a mathematical modeling experience using differential equations to model the volume of an intraocular gas bubble used by ophthalmologists to aid the healing of a surgically repaired region of the retina. We ask students to compare the traditional ``empirical'' model used by the ophthalmological profession with  more ``analytic'' models built on first principles and modeled with differential equations by examining the model predictions, after model parameters are estimated, using real patient data.  Finally, we introduce a new way of comparing models,  the Akaike Information Criterion, and apply it to our models.
\end{abstract}


\section*{NOTES FOR TEACHERS}




During some vitreoretinal surgeries the opthalmological surgeon will inject a gas, sometimes air, but often various mixtures of air and  sulfur hexafluoride ($SF_6$) or  perfluorocarbons ($C_nF_{2n+2}$), into the fluid within the vitreous cavity of the eye (the interior of the eyeball) to help the healing process. This creates an intraocular bubble  to keep pressure on the surgically treated area and provides a tamponade effect that causes functional closure of the retinal breaks, retinal detachment, macular holes, etc. See \cite{MichelsWilkinsonRice1990} for surgical procedure and details. 

\begin{quote}
Gas is now used quite frequently in many, if not most, operations; either approximately 16-20\% mixture of gases ($SF_6$ and air  or $C_3F_8$ and air) depending upon the case. Air, by itself, is very rarely used.\cite{Lambert2014}
\end{quote}

\noindent  By tamponade we mean the ``pathologic compression of a part.''\cite{MedDict2012} In effect, ophthalmologists seek to determine ``whether the intraocular bubble will tamponade the retinal break sufficiently long for a chorioretinal adhesion to form around the break''\cite{Thompson1989}.

We quote from the medical literature concerning the need to monitor the size of the bubble.
\begin{quote}
``Intravitreal  gas is commonly used in conjunction with vitrectomy and retinal reattachment to tamponade retinal breaks in the treatment of complicated retinal detachments. For a given retinal break, the size and duration of the intravitreal gas bubble determine whether effective tamponade of the retinal break is present sufficiently long enough to maintain apposition of the retina and retinal pigment epithelium. It becomes important to be able to accurately predict the kinetics of disappearance of intraocular gas bubbles to make an appropriate decision regarding the amount and type of gas to be used.''\cite[p.~609]{WongThompson1988} 
\end{quote}

Moreover, surgeons need to know what gas to use and how much gas to use. 

\begin{quote}
``Predictability of the kinetics of intraocular gases is important since the surgeon must choose the type of gas and the size of the intraocular gas bubble to achieve adequate duration of intraocular tamponade of the retinal break.''\cite[p.~691]{Thompson1989}
\end{quote}

The general shape of the eye is presumed to be a sphere of radius 1 cm and the gas bubble is known to have the following shape: a flat circular disk on the bottom (meniscus) with a domed top in contact with the eyeball. The bubble floats at the top of the eyeball and when the patient is standing upright this bubble is bounded below by the flat  disk meniscus, parallel to the ground. However, due to the optics of the eye, the bubble appears to the patient to be on the bottom of the eyeball with a flat top. The patient can see the exact shape of the bubble and its progress throughout the recovery period after surgery. 



Doctors monitor the size of the bubble in order to determine if the bubble is  suitable  to provide the tamponade effect. In Figure~\ref{fig:1-30-EyeWongThompsonEyebalP610} we see one approach to measuring the size of the bubble by recording the height  of the meniscus, i.e. the flat base of the bubble from the base of the eyeball, and then computing the height at time $t$ (in days) as a percentage of the vertical diameter of the eye. This is easily converted to centimeters using the assumption that the typical eyeball is 2 cm in diameter and so $h = h(t)$, the actual height of the meniscus in cm, is simply one minus the percentage  observed times 2 cm, i.e.

\begin{equation}\label{eq:height}
h = \left(1 - \frac{{\rm Percent}}{100}\right)*2\, . 
\end{equation}

\begin{figure}[h]
\centering

\includegraphics[width=7cm]{1-30-EyeWongThompsonEyebalP610.eps}



\parbox{12 true cm}{\caption{\label{fig:1-30-EyeWongThompsonEyebalP610}  ``The meniscus height of an intraocular gas bubble is estimated clinically as a percentage of the vertical diameter of the eye seen through the dilated pupil with the plane of the cornea perpendicular to the ground.''\cite[p.609]{WongThompson1988}}}
\vspace*{-.3 true cm}
\end{figure} 
           

Usually, the surgeon injects sufficient gas to occupy all of the spherical vitreous volume of the eyeball.  Thus we shall assume the bubble occupies  the entire volume of the eyeball (vitreous cavity) initially, i.e. $h(0) = 0$. 

Over a period of time the gas diffuses out of the eyeball and the full vitreous fluid level returns to the patient's vitreous cavity. A reasonable amount of time for such a bubble to go to a negligible level (say less than 0.1\% of the volume of the eyeball) is 60 days.   Interestingly, in patients with diabetes or other reasons that do not allow the intraocular fluid to exit normally, the bubble lasts much longer.\cite{Lambert2014} 



\subsection*{Patient Data}

In personal correspondence with Dr.~John T. Thompson \cite{Thompson1996} we received a  set  of human data which was used  in the study found in \cite{Thompson1989}. We have placed the complete data set in a supplemental file to this paper  under the name 1-030-EyeData.xls.  This file contains  an explanation of the exact nature of the data file, as well as the data itself. We selected several patient data sets for our work here and offer them in Table~\ref{table:PatientData}. 


\subsection*{Initial Model and Rationale}


The medical literature suggests that  the volume of the intraocular  gas bubble in the vitreous cavity, $V = V(t)$ in cm$^3$ at time $t$ in days can be  modeled by the following exponential decay model

\begin{equation}
\label{eq:initial}
\frac{dV}{dt} = - k V \, , \qquad  V(0) = \frac{4}{3}\pi (1^3)\, .
\end{equation}
NB:  The radius of the typical eye is presumed to be 1 cm.


Indeed, in calculating the absorption of the  bubble there are two models studied:  ``(1) natural logarithm of bubble volume vs.~time, and (2) surface area of contact of the intraocular bubble vs.~time. \ldots With model 1, the absorption of nonexpansile intraocular gas bubbles has been reported to approximate a first-order exponential equation with respect to intraocular gas bubble volume vs time \cite{WongThompson1988,LincoffEtAl1984,AbramsEtAl1974}.''\cite{Thompson1989} From a summary plot  the author observes ``Compared with model 2, model 1 shows a more rapid initial decline in absolute volume but a slower decline in absolute volume when the intraocular bubble is small.''\cite[p. 689]{Thompson1989} There is support for the exponential decay model in practice through these observations due to its ease of use and the use of the common notion of half-life for the height of the bubble.


We consider data from several patients, identifying them by patient number as found in  \cite{Thompson1996} for  analysis. In Table~\ref{table:PatientData} we offer the complete data set for these patients from \cite{Thompson1996}.


\begin{table}[h]\small
\begin{center}
\begin{tabular}{|c||c|c|c|c|c|r|l|l|l|l|} \hline \hline
126	&I	&C	&0.1	&6	&25	&1	&0.90	&4.5850	&1.5228	&0.2\\  \hline				
126	&I	&C	&0.1	&6	&25	&2	&0.90	&4.5850	&1.5228	&0.2\\  \hline					
126	&I	&C	&0.1	&6	&25	&3	&0.73	&3.8377	&1.3449	&0.54\\  \hline					
126	&I	&C	&0.1	&6	&25	&6	&0.45	&2.0266	&0.7064	&1.1\\  \hline					
126	&I	&C	&0.1	&6	&25	&13	&0.33	&1.2460	&0.2199	&1.34\\  \hline					
126	&I	&C	&0.1	&6	&25	&16	&0.25	&0.7830	&-0.2446	&1.5\\  \hline					
126	&I	&C	&0.1	&6	&25	&22	&0.10	&0.1460	&-1.9241	&1.8\\  \hline					
126	&I	&C	&0.1	&6	&25	&34	&0.02	&0.0062	&-5.0832	&1.96\\  \hline	\hline
179	&P	&C	&0.2	&12	&2	&1	&0.85	&4.8872	&1.5866	&0.3\\  \hline
179	&P	&C	&0.2	&12	&2	&2	&0.80	&4.6630	&1.5397	&0.4\\  \hline
179	&P	&C	&0.2	&12	&2	&3	&0.75	&4.3912	&1.4796	&0.5\\  \hline
179	&P	&C	&0.2	&12	&2	&6	&0.70	&4.0791	&1.4059	&0.6\\  \hline
179	&P	&C	&0.2	&12	&2	&10	&0.70	&4.0791	&1.4059	&0.6\\  \hline
179	&P	&C	&0.2	&12	&2	&20	&0.45	&2.2137	&0.7947	&1.1\\  \hline
179	&P	&C	&0.2	&12	&2	&36	&0.28	&0.9961	&-0.0039	&1.44\\  \hline
179	&P	&C	&0.2	&12	&2	&42	&0.28	&0.9961	&-0.0039	&1.44\\  \hline	\hline	
186	&A	&C	&0.2	&13	&4	&1	&0.75	&3.2540	&1.1799	&0.5\\  \hline
186	&A	&C	&0.2	&13	&4	&2	&0.75	&3.2540	&1.1799	&0.5\\  \hline
186	&A	&C	&0.2	&13	&4	&3	&0.70	&3.0223	&1.1060	&0.6\\  \hline
186	&A	&C	&0.2	&13	&4	&6	&0.65	&2.7730	&1.0199	&0.7\\  \hline
186	&A	&C	&0.2	&13	&4	&13	&0.55	&2.2338	&0.8037	&0.9\\  \hline	\hline	
\end{tabular}
\end{center}

\vspace*{-.1 true cm}
\caption{\label{table:PatientData} Observational data for several  patients referred to by patient number  in which we see, for example,   Patient  126  has an Intraocular lens (I);  gas used for the bubble was Perfluoropropane gas ($C_3F_8$) (C);  concentration of the gas was 10\% (0.1);  patient was assigned series number (6);  patient was assigned a number in that series for the study (25);   times of observation in days were 1, 2, 3, 6, 13, 16, 22, and 34; the percentage (decimal) of meniscus height of the bubble (beginning with 0.90); the computed bubble volume (cm$^3$) from meniscus height (beginning with 4.5850);   natural logarithm of bubble volume (beginning with 1.5228); and  actual height of the bubble in cm (beginning with 0.2) (added by this author).   Patients 179 and 186 have other conditions denoted by P and A respectively. }
\vspace*{-.3 true cm}
\end{table}


In Figure~\ref{fig:1-30-EyePatientMeniscusHeightPlot} we  plot the meniscus height of the intraocular gas bubble for Patient 126 as estimated clinically as a percentage of the vertical diameter of the eye seen through the dilated pupil with the plane of the cornea perpendicular to the ground. We use (\ref{eq:height}) to convert from percent of maximum height of 2 cm to actual height in cm.

\begin{figure}[h]
\centering

\includegraphics[width=10cm]{1-30-EyePatientMeniscusHeightPlot.eps}



\parbox{12 true cm}{\caption{\label{fig:1-30-EyePatientMeniscusHeightPlot}  Plot of the meniscus height of the intraocular gas bubble for Patient 126.}}
\end{figure}

Figure~\ref{fig:1-30-EyePatientDataPlot} shows the computed volume of the  bubble for the observed data for Patient 126    as a function of time in days and would appear to motivate the immediate conjecture of an exponential decay model, i.e. the rate of change in volume per unit time in days is proportional to the volume itself.

\begin{figure}[h]
\centering

\includegraphics[width=10cm]{1-30-EyePatientDataPlot.eps}



\parbox{12 true cm}{\caption{\label{fig:1-30-EyePatientDataPlot}  Plot of the volume  of the intraocular gas bubble   for Patient 126.}}

\vspace*{-1 true cm}
\end{figure}


\subsection*{Activity 0} 

Explain  formula (\ref{eq:height}) and then show how one obtains  the volume of gas of the bubble depicted in Figure~\ref{fig:1-30-EyeWongThompsonEyebalP610} given its height.

\vspace*{.3 true cm}

Now, in Activities (1) - (5) below perform your analyses using Patient 126 data and then take one other patient from the file 1-030-EyeData.xls  and perform the same analyses.

\subsection*{Activity 1}

Either by examining some of the data in 1-030-EyeData.xls  or reasoning and some assumptions (state them) defend why  (\ref{eq:initial})  is a reasonable model.



\subsection*{Activity 2}  Pick several patient's data from the file 1-030-EyeData.xls  and model the data with (\ref{eq:initial}).  Determine the rate constant $k$ for several patients and compare them. Indeed, the decay rate $k>0$ is determined for patient data under various lens scenarios and the resulting model, so pick identical and different scenarios or conditions for the patient.

\subsection*{Activity 3}

Think hard about this physiological situation. The gas in the center of the bubble goes nowhere, but the gas along the wall of the vitreous region in the eyeball could escape. With that in mind offer up a different model than (\ref{eq:initial}). Here too, pick several patient's data from the file 1-030-EyeData.xls  and model the data with (\ref{eq:initial}).  Determine the parameters in your model for several patients and compare them. Indeed, the parameters may be determined for patient data under various scenarios and the resulting model, so pick identical and different scenarios or conditions for the patient.



\subsection*{Activity 4}
Let us refine our model in Activity 1-3 by considering ALL the different ways the gas can escape the vitreous region of the eyeball and reflect that in a final revised model. With this new model in mind offer up a different model. Here too, pick several patient's data from the file 1-30-EyeData.xls and model the data with this newest model.  Determine the parameters in your model for several patients and compare them. Indeed, the parameters may be determined for patient data under various scenarios and the resulting model, so pick identical and different scenarios or conditions for the patient.


\subsection*{Activity 5}

Compare results from all your models and attempt to defend why the opthalmological profession uses the exponential decay model with half-life information in its patient analyses.

\teacher{
\section*{NOTES FOR TEACHERS}

We offer complete analyses of the data for Patients 126, 179, and 186 in a Mathematica notebook, 1-030-T-Mma-EyeModel.nb  and a pdf form of that same file as well as analyses of several other patient data sets.

Students will quite quickly come up with a better model (\ref{eq:betterinitial}) for understanding what is going on inside the eye in the following model:

\begin{equation}
\label{eq:betterinitial}
\frac{dV}{dt} = - k S^{3/2} \, , \qquad  V(0) = \frac{4}{3}\pi (1^3) \, , \qquad  S(0) = 4 \pi (1^2).
\end{equation} 

This is because the gas dissipates through the surface area.  This is true in the case in which the bubble stays spherical throughout. However, the bubble does not stay spherical; indeed, it is a sphere with a polar cap sheared off the bottom of the spherical bubble and the cap of vitreous liquid gets bigger and bigger until the volume of the bubble itself is 0.  We seek to model this more complex geometry and then compare this new model with the  exponential decay model used in practice.  

\subsubsection*{Model building -- some options}

In our studies (both with students in Spring 1995 and on our own since that time) we posed a number of different models for the size of the bubble as a function of time, principally modeling for the height, $h(t)$,  of the bubble in cm from the meniscus of the bubble which is at the bottom of the bubble to the top of the eye and thence computing the volume of the bubble, $V(t)$, in ${\rm cm}^3$.


Basically the issues are what kind of model are we willing to work for and will our model serve our needs? In the case of the {\it empirical\/} model effort to offer up the exponential decay model for the volume of the bubble as a function of time we can understand this because of the shape of its plot and the reasonable assumption that the bigger the bubble the more the bubble decreases per unit time and the smaller the bubble the less the bubble decreases per unit time. Thus if $V(t)$ is the volume, in cm$^3$, of the bubble at time $t$, in days, then perhaps (\ref{exp:equation}) is an appropriate model. 
\begin{equation}\label{exp:equation}
\frac{dV(t)}{dt}= - k V(t)\, , \quad V(0) = V_0
\end{equation}

\noindent Moreover, such a model has an attendant ease in computation with the notion of half-life readily available, possessing a  relatively simple calculation  for the length of time to halve, and halve, and halve, etc. the size of the gas bubble at which time  the physician can be confident there is no danger.  This is when the gas bubble will be sufficiently small, say, for the patient to     on airplanes and the wound could be considered healed.    Even then, a very small bubble would cause a significant rise in pressure if the patient was on an airliner and experienced rapid decompression.\cite{DickertEtAl1986} 

Here $k$ is the rate constant with units $1/t$, in our case 1/day.  However, (\ref{exp:equation}) does not adequately explain how the disappearance of gas actually takes place and maybe a model like (\ref{equation:FirstOrderSurfaceArea}), 
where $S(t)$ is the surface area of contact of the bubble with the retinal surface  at time $t$, is more reasonable for it ``explains'' how the gas is disappearing -- through the retinal area in contact with the bubble.  

\begin{equation}\label{equation:FirstOrderSurfaceArea}
\frac{dV(t)}{dt}= - k S(t)\, , \quad V(0) = V_0\, , \quad S(0) = S_0\, .
\end{equation}

\noindent However, (\ref{equation:FirstOrderSurfaceArea}) means our rate ``constant'', $k$, does {\it not\/} possess our usual units of 1/day, but rather must have units cm/day.   If we were to change (\ref{equation:FirstOrderSurfaceArea}) slightly to 

\begin{equation}\label{equation:ThreeHalvesOrderSurfaceArea}
\frac{dV(t)}{dt}= - k S(t)^{3/2}\, , \quad V(0) = V_0\, , \quad S(0) = S_0\, , 
\end{equation}

\noindent then the units on our rate constant, $k$, are 1/day.  We used a similar model to predict the size of a piece of dry ice as it sublimated, assuming the rate of change in the volume (and hence mass when volume is multiplied by the mass density) was proportional to surface area raised to the $3/2$ power. This was reasonable, for  the dry ice would sublimate on the surface and our units needed to be correct.\cite{Winkel2011, Winkel2014}

 
Moving further, (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) admits to the possibility that the gas may leave the bubble {\it both\/} through the  contact surface area with the retina at the top of the bubble {\it and\/} through the fluid in the eyeball at the base of the bubble. 

\begin{equation}\label{equation:ThreeHalvesOrderMeniscusSurfaceArea}
\frac{dV}{dt}= - k_1 S(t)^{3/2} - k_2 M(t)^{3/2}\, , \quad V(0) = V_0\, , \quad S(0) = S_0\,  , \quad M(0) = M_0\,.
\end{equation}



As the  surface area  of the meniscus' disk at the base of the bubble, $M(t)$, has units ${\rm cm}^2$ then $M(t)^{3/2}$ would have units ${\rm cm}^3$. This would make  our rate constant, $k_2$, for the term $k_2 M(t)^{3/2}$ in our model (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) have units 1/day with which we are  comfortable.

This certainly makes the model more complicated, admitting two rate constants to be determined, not just one.  However,  (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) (where $M(t)$ is the surface area of the meniscus at the bottom of the bubble) would make more sense, for gas could escape through the meniscus surface area into the vitreous fluid which is increasing in the eyeball as well as through the  retinal surface area.  


 Thus we have articulated the two ways in which gas can escape, (1) the surface area of the retinal wall of the bubble and (2) the portion of the bubble  in contact with this meniscus's  surface area, not with the retinal wall, but rather with the vitreous fluid in the eyeball. Here, $k_1$ and $k_2$ are rate constants, both with the units 1/day. Thus, the general model in (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) just might make us most happy.  


Of course, we could just try a general power law model  (Model (3))  such as

\begin{equation}\label{equation:PowerLawSurfaceArea}
\frac{dV}{dt}= - k S(t)^{r}\, , \quad V(0) = V_0\, , \quad S(0) = S_0\,  .
\end{equation}

\noindent Here  (\ref{equation:PowerLawSurfaceArea}) would force our rate ``constant'' $k$ to have units ${\rm cm}^{3-r}$/day, another non-desirable situation.

Thus it would appear that model (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) has the inside track for an {\it analytic\/} model, i.e. one which makes a good attempt at explaining the reason (the mechanism) for the disappearance of the gas as well as possessing good units. 

 

Implicit in the ophthalmological literature is the fact that any diffusion of the gas across the meniscus of the growing vitreous fluid volume is negligible when compared to the absorption of the gas into the  surface area of contact on the retinal surface of the bubble. Because of surface tension at the interface of the vitreous fluid and the gas no gas will enter the liquid, for the overall pressure on the gas is insufficient to force it into the liquid. As every carbonated beverage drinker knows, gases seek to escape liquids rather then ``push'' their way in, as the carbon dioxide gas is forced to do in the production of carbonated beverages.  

However, in (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) we decided to allow the diffusion of the gas into the  retinal wall as well as into the vitreous fluid  through the meniscus surface area as well to see if the data supported this implicit assumption. We were then estimating two parameters $k_1$ and $k_2$ in model (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}).

\subsection*{Parameter estimation approaches}

 We could estimate the two parameters $k_1$ and $k_2$ through a process similar to that which we have used before \cite{Winkel2011a}. Namely, we create a reasoned grid in our $(k_1,k_2)$ parameter space and compute the sum of square errors between the timed observed data for the height of the meniscus ($h(t)$) and the predictions of   our model at these times using our numerically computed solution to (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) for each pair of   $k_1$ and $k_2$, i.e.

\begin{equation}\label{eq:SSE}
{\rm SSE}(k_1,k_2) = \sum_{i=1}^n (h_i - \hat{h_i}(t_i))^2 \, ,
\end{equation} 
where $n$ is the number of observations for a given patient, $h_i$ is the  observation of bubble height at time $t_i$, and $\hat{h_i}(t_i)$ is the model bubble height at time $t_i$.  





We then sort through all the computed sum of square errors for the minimum and use the corresponding values of $k_1$ and $k_2$ as our best parameters for our best model.  


\vspace*{.2 true cm}

Here we actually offer a method in which we use a gradient search algorithm on the $(k_1, k_2, {\rm SSE}(k_1,k_2))$ three dimensional parameter surface to determine the parameters $(k_1, k_2)$ which give the lowest ${\rm SSE}(k_1,k_2)$ and thence make our model fit the data best. 

We seek to determine the parameters $(k_1, k_2)$ which minimize the sum of square errors between our differential equation model and our data,  ${\rm SSE}(k_1,k_2)$. This sum of square errors, ${\rm SSE}$,  will be a function of the two parameters $(k_1, k_2)$, something like ${\rm SSE} = {\rm SSE}(k_1,k_2)$.  

We  give an initial estimate of our parameters $(k_1, k_2)$ and then compute nearby points in three dimensional parameter space $(k_1, k_2, {\rm SSE}(k_1,k_2))$ from which we will compute numerical estimates of the partial derivatives of our model's sum of square errors, ${\rm SSE} = {\rm SSE}(k_1,k_2)$,  as a function of the parameters $k_1$  and $k_2$ (using symmetric differences) in both $k_1$  and $k_2$ directions. 

For a given point in space on the surface ${\rm SSE} = {\rm SSE}(k_1,k_2)$ we then compute the gradient and take steps in the negative gradient which we know will decrease our ${\rm SSE}$ value most at that point. 

We then move a step in that direction and recompute the value of ${\rm SSE}$ at the new point, comparing it to the value of ${\rm SSE}$ at the previous point.  If it is less then we jump to the new point. If it is more then we cut our step size by a given factor  to recompute the gradient and take a new step to a closer point on our parameter surface  ${\rm SSE} = {\rm SSE}(k_1, k_2)$. 

We continue in this manner for a large number of iterations  and we plot our results in two parameter spaces, $(k_1, k_2)$ to see the convergence to our minimal ${\rm SSE}(k_1,k_2)$. 

\subsubsection*{Let the modeling begin}

We pursue the exponential decay model, but we also pursue our various models   (\ref{equation:ThreeHalvesOrderSurfaceArea}) - the rate of change in volume per unit time in days is proportional to $r=3/2$ power of the surface area of the bubble which is in contact with the retinal wall; and (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea})   -   the rate of change in volume per unit time in days is proportional to the sum of the rates at which  the gas  leaves the bubble {\it both\/} through the  contact surface area with the retina at the top of the bubble {\it and\/} through the fluid in the eyeball at the base of the bubble. 

We shall  use model (\ref{equation:ThreeHalvesOrderSurfaceArea}), to determine the rate constant parameter $k$ and compare the results with the exponential decay model approach taken in the papers \cite{Thompson1989,WongThompson1988}.  This we shall do using  patient data provided to us \cite{Thompson1996}.  Specifically, here we use data from Patient 126 as presented in Table~\ref{table:PatientData}.  
 
As a more academic, and challenging exercise, we   also address the broader model in (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}). We  could examine (\ref{equation:FirstOrderSurfaceArea}), for students always like to do (and actually propose) this simpler calculus before complicating their analysis with something like 
(\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}). Also, we  could consider 
(\ref{equation:PowerLawSurfaceArea}) in the simple case where $r = 2$ because, again, students go there in their own analysis as it appears easier than more complicated models involving fractional powers (\ref{equation:ThreeHalvesOrderSurfaceArea}) or  another outlet for the gas found in 
(\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}). Finally, we consider the most general model (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}) and estimate its parameters while comparing its predictive value to that of our previous models. 
  

We note in Figure~\ref{fig:1-30-EyePatientThreeHalvesExpPlots} that there is little difference in the plots of the   models, Model (1), i.e. Exponential Decay Model (thick curve) and Model (2) $r = 3/2$ (thin curve) for Patient 126 data. Moreover, from Table~\ref{table:Summary} the Sum of Square Errors (SSE) between the model and the data for the Model (1)  is  0.267657 while with Model (2) it is 0.283921. Even the numbers do not differentiate much. However, Model (3), our general power law,  does differ with  Model (1) and Model (2), both in how it appears to fit the data better (see Figure~\ref{fig:1-30-EyePatientThreeHalvesExpPlots}) and it is considerably smaller SSE, i.e., 0.101588.


\begin{figure}[h]
\centering

\includegraphics[width=10cm]{1-30-EyePatientThreeHalvesExpPlots.eps}



\parbox{12.8 true cm}{\caption{\label{fig:1-30-EyePatientThreeHalvesExpPlots}  Summary plot of best fit models for  Patient 126 data, Model (1)  (thick curve), Model (2) (thin curve), and Model (3) (dashed curve).  }}  

\vspace*{-.4 true cm} 
\end{figure}

\section*{Comparing Models}


Thus we see that we  have several models for the volume of the intraocular bubble of gas in   the eye. The first, (\ref{equation:FirstOrderSurfaceArea}), has one parameter $k$;   the second, (\ref{equation:ThreeHalvesOrderSurfaceArea}), has one parameter, $k$; and the third (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}), has two parameters $k_1$ and $k_2$. The first has a sum of square errors of 0.267357;   the second has a sum of square errors of 0.283921; and the third has a sum of square errors of 0.101588. 

Some texts refers to the sum of square errors (SEE) as the {\it residual sum of errors\/} (RSE). Now certainly if we increase the number of parameters we can get a better fitting model in terms of reducing the SSE. Indeed, if we took an 8$^{\rm th}$ degree polynomial we could fit it to go through all the 8 data points and our SSE would be 0. However, this model would not be valuable in any predictive sense. Moreover, it would not explain anything as to the actions of the underlying model. So there is a trade-off, namely,  if we increase the number of parameters in our model we can decrease our SSE. We can ask which model  is best, in general. Namely which model keeps both the SSE and the number of parameters low.

In 1974 a formula  was published  \cite{Akaike1974} which would allow us to use both SSE and number of parameters at the same time to compare models. The statistic is called the {\it Akaike Information Criterion\/} or AIC.  A reasonable formulation (\ref{eq:AIC}) for our purpose is given in  \cite[p.~138]{Ledder2013}:


\begin{equation}\label{eq:AIC}
{\rm AIC}(n, k, {\rm SSE}) = n \ln\left(\frac{{\rm SSE}}{n}\right) + 2 (k+1) \,,
\end{equation}
where $n$ is the number of observations, $k$ is the number of parameters to be estimated, and SSE is the least sum of square errors.

AIC is used to compare models and the model with the smallest AIC is considered best of the set if the difference between that model's AIC and all other AIC's is more than 2. Otherwise the AIC cannot tell us which model is best.  AIC does not tell us if the model is good, only if it is the best of the ones being considered.

Smaller (even negative) values of AIC represent a higher level of statistical support for the corresponding model. The smallest AIC indicates the best model - in terms of balancing ACCURACY (minimum sum of square errors) and COMPLEXITY (number of parameters).  A difference of 2 in the AIC's is needed to suggest a significant difference. 
 
If all AIC's are close to each other then Model (1), the exponential fit, is just as good as the more exhaustive Model (3), using the surface area of the bubble and the surface area of the membrane of the meniscus. 

We apply  this AIC analysis to our situation.


\section*{SUMMARY OF ANALYSES}

We offer a summary (see Table~\ref{table:Summary}) of three patients, Patients 126, 179, and 186. While in all three cases Model (3) offers the lowest SSE    only  Model (3) in the cases of Patients 126 and 186 possess an AIC value  sufficiently below (more than 2) the other two models to say Model (3) is the best of the three models.  However, in the case of Patient 179  all the models  seem to do a good job and the graphics in all cases seem to be comparable to those obtain for Patient 126 depicted in Figure~\ref{fig:1-30-EyePatientThreeHalvesExpPlots}.



\begin{table}[h]


\begin{center}\small 
\begin{tabular}{||c|c|c|c|c|r||} \hline
Patient  &Model & \multicolumn{1}{|c|}{$k$}& \multicolumn{1}{|c|}{SSE}&\multicolumn{1}{|c|}{AIC }\\ 
Number & \multicolumn{1}{|c|}{Number} &  \multicolumn{1}{|c|}{Constant(s)}& &\multicolumn{1}{|c|}{Value}\\ \hline \hline 
126	&1&\multicolumn{1}{|c|}{0.10896}&0.267357&-23.1889\\ \hline 
126	&2&\multicolumn{1}{|c|}{0.01365}&0.283921&-22.7080\\ \hline  
126	&3&\multicolumn{1}{|c|}{0.00992,\,\, 0.017031}&0.101588&-28.9301\\ \hline \hline 
179	&1&\multicolumn{1}{|c|}{0.04082}&0.048482&-36.8480\\ \hline 
179	&2&\multicolumn{1}{|c|}{0.00510}&0.046901&-37.1132\\ \hline  
179	&3&\multicolumn{1}{|c|}{0.00330,\,\, 0.005824}&0.034038&-37.6777\\ \hline \hline 
186	&1&\multicolumn{1}{|c|}{0.05822}& 0.081213&-16.6006\\ \hline 
186	&2&\multicolumn{1}{|c|}{0.006571}&             0.065483&-17.6770\\ \hline  
186	&3&\multicolumn{1}{|c|}{0.00330,\,\, 0.005169}&0.002198&-32.6480\\ \hline \hline
\end{tabular}




\parbox{11.9 true cm}{\caption{\label{table:Summary} Summary table of data for Patients 126, 179, and 186 in which we display Model Number, rate constant(s) $k$'s, Sum of Square Errors (SSE), and AIC Value.}}
\end{center}
\end{table}

Now, from our computed AIC values for the three models for Patient 126 as given in Table~\ref{table:Summary} it is clear that Model (3) is the best model. For Patient 126 we have an AIC for our first model, (\ref{equation:FirstOrderSurfaceArea}), of -23.1889; an AIC for our second model, (\ref{equation:ThreeHalvesOrderSurfaceArea}), of -22.708; and an AIC for our third model, (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}), of -28.9301. Note, the AIC for Model (3) is at least 2 lower than that for Models (1) and (2).  This means our best model, according to the AIC,  is the more complex model  (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea})  using the surface area of the bubble and the surface area of the membrane of the meniscus.  However, recall that  in the ophthalmological literature  the fact that any diffusion of the gas across the meniscus of the growing vitreous fluid volume is negligible when compared to the absorption of the gas into the  surface area of contact on the retinal surface of the bubble, i.e., $k_2 = 0$. Thus, this would have us not really  accept  Model (3) as totally reasonable for $k_2 = 0.017031 \neq 0$ and we just might be willing to settle for the profession's choice, namely,  the exponential decay model, Model (1), after all.

Indeed,  what can we say after posing alternate models to the profession's immediate use of Model (1)? Certainly, such a model is quite natural given an examination of the plot of the data. See Figure~\ref{fig:1-30-EyePatientDataPlot} and  note that all patient data looks quite similar in nature. The intuitive notion that the bigger the bubble is the more volume it will lose per unit time simply because of the exposure of the bubble to more {\it surface area\/} of the retinal wall seems to hold up. However, it is exactly those italicized words in the previous sentence, namely {\it surface area\/}, which makes us still prefer Model (2)  as being  better at explaining the data. Additionally, because of the correctness of the units involved in Model (2), i.e. $r=3/2$,  we believe this model is the best at explaining the phenomena, while at times not offering the smallest  sum of square errors or a significant difference in AIC. Indeed, the ability to determine half-life with the exponential decay model and the importance of quickly assessing whether or not the bubble will last long enough to do its healing purpose makes the exponential decay model, Model (1), the preferable one for the profession. Moreover, it works quite well as a bubble volume predictor and measure of safety and healing power of the whole tamponade effect in post-surgery retinal situations.



 






\section*{CONCLUSION}


In the winter of 1994-95 I experienced a torn retina and I was successfully treated and operated on by Dr.~Hector Ibanez of Indianapolis, Indiana. As a result of this experience I became interested in one aspect of the surgical procedure.  I am grateful to Dr.~Ibanez for providing me with copies of medical text materials \cite{MichelsWilkinsonRice1990} so I could better understand the surgery approaches, the geometry and physiological processes of the eye, and the experience I went through in successful recovery.

In all patient data we examined we noted that with all three models offered, namely (\ref{equation:FirstOrderSurfaceArea}), (\ref{equation:ThreeHalvesOrderSurfaceArea}), and (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}), approximately the same best fitting model was obtained and so we would concur with the ophthalmological profession that the  use of the exponential decay model is warranted. For there is no further information to be gained from the two more complex models, namely, (\ref{equation:ThreeHalvesOrderSurfaceArea}) and (\ref{equation:ThreeHalvesOrderMeniscusSurfaceArea}).

In this scenario we took a personal health experience and found a mathematical modeling opportunity for our students in the situation. We shared  with our students who built and rejected models (we were part of the data!), asked questions, solved equations, and validated their models. We learned of the simple exponential decay model used in ophthalmological practice and thought that a more reasoned model, based on a deeper understanding of  the physics of the situation in the eye, would serve better. Our studies found that this was not always the case, at least based on comparisons of various models using a least sum of square errors criterion or an Akaike Information Criterion comparison for the patient data sets we used. This venture into real world modeling proved good for our students and for us, perhaps that ``us'' will include you as you find models in and around your own life. 





\section*{ACKOWLEDGEMENT}

Of course, we are grateful to our surgeon, Dr.~Hector Ibanez, Indianapolis IN, for the successful surgery and for introducing us to the ophthalmological literature through a ``pile'' of copied material from journals and articles. Indeed, we are very grateful for his tact  in our consultations  for an eye problem as he would say that unless the gas has fully dissipated one could have difficulties in air travel in the event of a loss of cabin pressure.  We are also grateful to Dr. John T. Thompson of The Retina Institute of Maryland for sending the data set pertaining to the human data published in  \cite{Thompson1989}.

We are also grateful to Dr.~Michael Lambert, an ophthalmological surgeon, whose experiences include being an officer in the United States Air Force  and serving as a test pilot and physics faculty member at the United States Air Force Academy, for his many suggestions on the presentation and technical details.

}

\begin{thebibliography}{1}

\bibitem{AbramsEtAl1974} Abrams, G. W. , H. F. Edelhauser, T. M. Aaberg, and L. H. Hamilton.  1974.  Dynamics of intravetireal sulfur hexafluoride gas.  {\it Investigative Ophthalmology \& Visual Science\/}. 13: 863-868.
 
\bibitem{Akaike1974} Akaike, H. 1974.  A new look at the statistical model identification. {\it  IEEE Transactions on Automatic Control\/}.  19(6): 716-723.


\bibitem{DickertEtAl1986}  Dieckert, J. P., P. J. O'Connor,   D. S. Shacklett, T. J. Tredici, H. M. Lambert, J. Fanton, and J. O.  Sipperley.  1986. Air Travel and Intraocular Gas.  {\it Ophthalmology\/}.  93(5): 642-645. 

\bibitem{MedDict2012} {\it  Free Dictionary by Farlex\/}.  \href{http://medical-dictionary.thefreedictionary.com/tamponade}{http://medical-dictionary.thefreedictionary.com/tamponade}. Accessed 12 September 2014.

\bibitem{Lambert2014} Lambert, M. Personal correspondence. 16 March 2014.

\bibitem{Ledder2013} Ledder, G. 2013. {\it Mathematics for the Life Sciences\/}. New York: Springer. 

\bibitem{LincoffEtAl1984} Lincoff, H. A., J. M. Maisel, and A. Lincoff. 1984. Intravitreal disappearance rates of four perfluorocarbon gases.  {\it Archive of ophthalmology\/}. 102: 928-929.


\bibitem{MichelsWilkinsonRice1990} Michels, R. G., C. P. Wilkinson, and T. A. Rice. 1990. {\it Retinal detachment\/}. St. Louis MO:  The C. V. Mosby Company.


\bibitem{Thompson1989} Thompson, J. T. 1989. Kinetics of Intraocular Gases:  Disappearance of Air, Sulfur Hexafluoride, and Perfluoropropane After Pars Plana Vitrectomy. {\it Archive of ophthalmology\/}. 107: 687-691.

\bibitem{Thompson1996} Thompson, J. T. 1996. Personal correspondence. Letter dated 9 March 1996.


\teacher{
\bibitem{Winkel2011a}Winkel, B. J. 2011. Parameter Estimates in Differential Equation Models for  Population Growth. {\it PRIMUS\/}.  21(2): 101-129.

\bibitem{Winkel2011} Winkel, B. J. 2012. Modeling Sublimation of Carbon Dioxide.  {\it International Journal of Mathematical Education in Science and Technology\/}. 43(8): 1077-1085.  


\bibitem{Winkel2014}  Winkel, B. J. 2014. 1-012-T-Sublimation of Carbon Dioxide.  
 

\bibitem{Winkel2014a}  Winkel, B. J. 2014.  Modeling Intraocular Gas Bubbles in Retinal Surgery Patients -- A Differential Equation Approach. Submitted for publication.

}

\bibitem{WongThompson1988}
Wong, R. F. and J. T. Thompson. 1988.  Prediction of the Kinetics of Disappearance of Sulfur Hexafluoride and Perfluoropropane Intraocular Gas Bubbles.  {\it Ophthalmology\/}. 95(5): 609-613.


\end{thebibliography}




\end{document}
