Establishing Similarity of Modeled and Experimental PK/PD Hysteretic Loops using Pseudo Time Series Analysis and Dynamic Time Warping Ron Smith, Florida Southwestern State College, Fort Myers, Fl.
Abstract An iterative Bezier curve cubic spline algorithm is used to approximate a hysteretic loop for a pharmacodynamics(PD) dose/response curve. The accuracy of the resulting vector-valued function is tested using graphs produced through PROC SGPLOT. The centroid of each curve is found through PROC IML. Euclidean distances from each centroid to points on the respective polygonal curves (experimental and model) produce several distance vectors for each curve. A sequence of these distance vectors is plotted as a Pseudo Time Series for both the original experimental and the adapted model curves. The graphical similarity of these curves is examined using the Dynamic Time Warping (DTW) algorithm in PROC SIMILARITY. The information produced from the DTW process gives a numerical measure of how well the model follows the experimental curve. The parametric equations from this a posteriori model can form a composition with a time-based function to produce the final PK/PD effect versus time equations.
Introduction Hysteresis loops occur from the time lag of a response (effect) behind the instigator of the “effect”. In dose response curves the two typical reasons for the lag are (1) a dearth of drug receptor sites and (2) slowed drug-receptor site interaction. The aforementioned reasons produce counter-clockwise hysteresis loop, in which plasma drug concentration precedes then lags the observed effect with time. An alternative hysteretic curve process would include any situation in which the value of one variable depends on whether the other variable is increasing or decreasing. A clockwise hysteresis path occurs with the development of tolerance to a drug. Hysteresis loops are graphical traces occasionally encountered in the analysis of pharmacokinetic and pharmacodynamic relationships. Hysteretic curves often indicate more complex drug-receptor interaction and/or the presence of confounding drug metabolites. They are typically relations and not functions, forcing a different modeling process that might involve systems of equations.
Figure 1. Clockwise and Counter-clockwise Hysteresis loops.
Figure 1. (a) counterclockwise movement of a dose response curve where the drug and receptor interactions are limited. (b) clockwise movement where drug tolerance has increased. Pharmacokinetics (PK) measure effect versus time, while pharmacodynamics (PD) measure concentration versus effect. Ultimately a PK/PD curve will measure effect versus time (Fig. 2).
Pharmacokinetics (PK) [concentration] vs. time
Pharmacodynamics (PD) [concentration] vs. effect
PK/PD Effect vs. Time
Figure 2. The com bination of the dim ensions of the pharm acokinetics and pharm acodynam ics in pharm aceutical dose -response studies.
The present paper will describe an a posteriori model of the effect versus drug concentration hysteresis loop using a Bezier type curve adaptation that will produce two parametric equations of the form x( ) and y( ), with as a parametric placeholder for a function of time. That is, we have a composition where = f(t). Figure 3. shows how a time component may be displayed along with the planar hysteresis loop in an experimental setting. Hour 5 Hour 4
Hour 7 Hour 8
Hour 10 Hour 2
Hour 11 Hour 0
Drug Concentration Figure 3. Hourly assay points superim posed on a typical experim ental hysteresis plot. The dots represent tim e values w here blood assays are taken along the concentration effect curve. Note: Tim e is often not in linear relation w ith the unit m ovem ent along the curve. Tim e and distance along the curve m ay be logarithm ically or exponentially related.
There are many researchers who have developed a priori models for hysteretic systems. The equations for these systems are very briefly described in Appendix B.
The outline of the a posteriori process of the current paper follows: 1. A cubic Bezier spline curve is produced using an iterative process that determines the distances of 1/3 and 2/3 the measure along the complete curve path. Control points for the model curve are derived from points on the experimental curve that might be obtained through pharmaceutical research. 2. Once the centroid of each of the two curves is found, the Euclidean distance from the center of mass to each of the sequential points on the curve is found. The included Bezier Curve analysis insures that the starting and ending points (identified as anchor points) for the two curves are automatically the same. 3. A pseudo time series for each curve is determined using the Euclidean distance plotted against an integer n-value starting at 1 and ending with the nth point on each polygon. Here, the (n+1)st point is equivalent to the 1st point. 4. The pseudo time series are evaluated by PROC ARIMA and checked for ARIMA model and stationarity. 5. The pseudo time series values are examined in PROC SIMILARITY using Dynamic Time Warping to give a measure of the similarity of the two curves. A small difference found through DTW is taken as a confirmation that the curves are similar and that the adaptive curve might be used to model the experimental curve.
Modeling the Dose-Response Hysteretic Loop The shape of the hysteretic loop in the PK/PD dose-response curve is mindful of a cubic spline curve, with the cubic Bezier curve being a special form of the spline curve. A spline is a series of connected polygonal segments. The segments can be linear, quadratic, cubic, or even higher order polynomials. As we produce our model curve, we derive a set of parametric-form vector-valued functions using a set of four control points (Po, P1, P2, P3) that produce the cubic Bezier polynomial. The cubic Bezier curve is given by is found from the binomial expansion of Equation 1. Equation 1: Cubic Bezier (3, ) =
P ( ) (1 ) k
using the parameter, , and its complement (1- ) since [0,1] . The control points help form the curve, with Po, P3 as anchor points which are on the Bezier curve (black) and P1, P2 are control points not on the curve. The solid back dots are auxiliary control points with point A on the Bezier curve.
P3 Figure 4.The four control points of the cubic Bezier spline w ith interior tangent at point A.
Since points Po, P3 are already known, we need find points P1, P2 to produce the equation for the model curve. Po, P3 will be the same for the experimental and the adapted model curve. In order to get the best
fit for the model curve, we must obtain two points from the original curve for the dose-response data. Measuring (counter-clockwise) along the curve, it is common to select a point one-third of the distance from the anchor point, P o, to the terminal anchor point, P 3 and another point at two-thirds of the distance. In terms of a matrix equation, Equation 1 becomes (now identified as Equation 2): Equation 2:
B( ) T C P where B contains the two parametric equations, x( ) and y( ), that will trace the curve on the twodimensional plane. T is the power matrix, C is the coefficient matrix and P is the matrix of control points. In the expanded form, the matrices become: Equation 3:
x( ) 2 y ( ) 1
1 0 0 3 3 0 3 3 6 3 1 3 3
0 Po ( xo , yo ) 0 P1 ( x1 , y1 ) 0 P2 ( x2 , y2 ) 1 P3 ( x3 , y3 )
We need four data values from the experimental original curve and determine control point information that will interpolate our data for a model. We can state that we have four values Q0, Q1, Q2 and Q3, from our experimental curve and that we want to determine control values so that our Bezier function will pass through these data points. To force interpolation, it is only necessary to set the definition of the control point values as follows: Equation 4
27 0 0 0 1 8 12 6 1 Q P 27 1 6 12 8 0 0 0 27 The above formulas were derived by evaluating the formal representation of the cubic Bezier function at path distance = 0, 1/3, 2/3 and 1. The four equations produce the x-values and y-values of the control points we seek. The Bezier function P(n) determined by these control values Q ( ) will have the interpolation property we desired, namely, that: Equation 5 Q (0) = P 0 Q (1/3) = P 1 Q (2/3) = P 2 Q (1) = P 3 The relationship we derived above may be written as that found in Equation 6. This relationship is found through the inverse of the coefficient matrix:
6 0 0 0 1 5 18 9 2 P Q 6 2 9 18 5 0 0 0 6
The fact that each row contains positive entries that add up to 1.0 is an illustration of the fact that the data values are convex combinations of the control values. The value from the experimentally derived curve are as follows for twenty -one points identified. point 1 2 3 4 5 6 7 8 9 10 11
x 50 60 75 112 120 160 180 215 220 240 255
y 58 60 70 75 77 80 81 82 81 80 70
point 12 13 14 15 16 17 18 19 20 21
x 290 300 350 340 300 250 229 200 100 0
y 63 49 39 33 30 21 20 17 10 0
Table 1 Points Attained from Experim ental Curve . Num bered points (in shaded orange colum ns) w ith x and y values.
An in-depth examination of the curve’s arc gives Q 0 (50,58); Q1 (distance = 1/3) (220,81); Q2 (distance = 2/3) (350,39); Q3 (0,0). Multiplication of the four points by the coefficient matrix gives: Equation 7
xo x P 1 x2 x3
yo 6 0 0 0 50 y1 1 5 18 9 2 220 y2 6 2 9 18 5 350 y3 0 0 0 6 0
58 50 81 109.3 39 729.2 0 0
58 131.7 19.3 0
Which gives the two parametric equations: Equation 8
x( ) 50 174 1689 2 1913 3 y ( ) 58 219 555 2 278 3 Rewriting these parametric equations with the previously mentioned f(t) = , we can now plot the location on the curve with respect to the parameter time using
x( f (t )) 50 174( f (t )) 1689( f (t )) 2 1913( f (t )) 3 y ( f (t )) 58 219( f (t )) 555( f (t )) 2 278( f (t )) 3 These equations produce the graph shown below for the adaptive model Bezier curve. The numerical values associated with the dots on the curve are - values.
Figure 5. The Adapted Model Curve. Show ing param etric tau-values at 0.05 values along the curve. See Appendix C for SAS code.
From these parametric equations, we can add the modeled-curve graph to that of the original experimentally obtained curve below:
Figure 6. Graph of the Original Curve and the Adapted Model Curve. Com m on plot of the experim ental curve (red) and the m odeled curve (blue).
Point 1 2 3 4 5 6 7 8 9 10 11
x 50 62 75 107 137 169 202 236 267 296 320
y Point x y 58 12 338 57 68 13 349 50 74 14 351 42 79 15 343 35 81 16 323 27 82 17 291 20 81 18 243 14 79 19 180 8 75 20 99 3 69 21 0 0 63
Table 2. Num bered points (in shaded orange colum ns) w ith x and y values (the closed polygon actually starts at (0,0) and returns to (0,0)).
Produce a Centroid for Each Curve using PROC IML. The area of a simple closed polygon with vertices (x0, y0), (x1, y1), ..., (xn-1, yn-1), is given by Equation 10
1 n 1 ( xi yi1 xi1 yi ) 2 i 0
In the formula, the vertex (xn, yn) is identified with the first vertex, (x0, y0). The centroid of the polygon occurs at the point (cx, cy), where Equation 11
1 n 1 ( xi xi1 )( xi yi1 xi1 yi ) 6 A i 0 1 n 1 cy ( yi yi1 )( xi yi1 xi1 yi ) 6 A i 0
The above equations are employed in finding the “center of mass” (centroid) of closed polygons formed from the joining of the ends of the two curves under study (the experimental drug curve and the model cubic Bezier curve). Rick Wicklin, of the SAS Institute, has produced (in his blog, The Do Loop) a very fine SAS procedure in PROC IML to find the centroid (barycenter) of a simple closed polygon. This procedure (see appendix) was employed to find a barycentric point for each of the experimental curve and the adaptive model curve. The results of the aforementioned procedure are presented in the following tables and graphs. Experimental Drug Polygon Centroid x Centroid y
Table 3. Calculated x and y values of the centroid based upon the program furnished by Dr. Rick Wicklin in The DO Loop.
Dashed line completes
Figure 7. Plot of the experim ental curve show ing the calculated centroid position.
Adaptive Model Polygon Centroid x Centroid y 172.895 41.5548 Table 4.Calculated x and y values of the centroid based upon the program furnished by Dr. Rick Wicklin in The DO Loop.
centroid Dashed line completes polygon.
Figure 8. Plot of the m odeled curve showing the calculated centroid position.
Compared to the absolute distances found in the experimental and adapted model curves , the centroids are quite close in location. While the goal of the procedure is to test the similarity of the two curves, it is of note that the Euclidean distance between the centroids (although not exactly coincident) is proximal (approximately 9.1 units apart).
Finding the Euclidean distance from the center of mass to each of the sequential points on the curve. Using the distance formula, d
x2 x1 y2 y1 2
, the Euclidean distance from the centroid to each
point in the data set that lies on the curve is found. The Euclidean distance is the distance metric commonly used in simpler shape analyses. The results are in the two tables that follow: Point
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
50 60 75 112 120 160 180 215 220 240 255 290 300 350 340 300 250 229 200 100 0
58 60 70 75 77 80 81 82 81 80 70 63 49 39 33 30 21 20 17 10 0
Euclid Distance from Centroid 112.603 103.016 90.733 59.081 53.794 37.162 42.343 66.215 69.759 86.763 97.277 129.985 138.553 188.457 178.690 139.015 91.083 71.192 46.317 69.807 167.173
Table 5. Euclidean distances of points on the experim ental curve is found using the distance form ula.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
50 62 75 107 137 169 202 236 267 296 320 338 349 351 343 323 291 243 180 99 0
58 68 74 79 81 82 81 79 75 69 63 57 50 42 35 27 20 14 8 3 0
Euclid Distance from Centroid 123.991 114.005 103.132 75.792 53.333 40.632 49.020 73.378 99.871 126.127 148.660 165.825 176.307 178.105 170.231 150.809 120.055 75.325 34.299 83.349 177.819
Table 6. Euclidean distances of points on the m odeled curve is found using the distance form ula.
A Pseudo Time Series is formed.
A pseudo time series is determined using the Euclidean distance plotted against number of the nth sequential point on the curve. Here, the (n+1)st point is equivalent to the 1st point. Wei (et al.) demonstrates in the following figure how a plot of the Euclidean distance against the point designation produces a Pseudo Time Series.
Figure 9. Example used to demonstrate mapping of the centroid and Euclidean distance to convert a closed simple curve to a Pseudo Time Series. From Wei (2008).
Figure 10. The colored coded arrow s in the figure above point to the place in the Pseudo Tim e Series that relates to the Euclidean distance from that exam ple point to the centroid.
The process used in producing the time series for the experimental drug curve is also used to produce a Pseudo Time Series for the adaptive model curve
Figure 6. Plot of the Euclidean Distance versus the index num ber of the Pseudo Tim e Series.
The Pseudo Time Series Tested for Stationarity using PROC ARIMA The Euclidean distances for the Pseudo Time Series were obtained from a closed polygon and that process would necessarily indicate a stationary time series. However, the stationarity of each Pseudo Time Series was further tested through an ARMA model using PROC ARIMA. The following tables display the mean and standard deviation of the respective pseudo time series. Mean and Std. Dev. Model Mean of Working Series 111.4317 Standard Deviation 46.86071 Number of Observations 21 Table 7. Mean and St. Dev. of the Model Pseudo Tim e Series Found through PROC A RIMA
Mean and Std. Dev. Experimental Mean of Working Series 97.0961 Standard Deviation 43.9356 Number of Observations 21 Table 8. Mean and St. Dev. of the Experim ental Pseudo Tim e Series Found through PROC ARIMA
Example code for the model follows: proc arima data=adaptive; identify var=euc_dist(1) stationarity=(ADF=(0,1,2,3)) outcov=adf; label euc_dist = "Adaptive Model Curve"; run; The models for both times series shows that the MA (1) model is appropriate. This model assures us that the two Pseudo Time Series are stationary (since these times series are necessarily finite). If the respective time series are not stationary (i.e. with a constant mean), we would not be able to compare them properly.
SCAN (from PROC ARIMA) Experim ental Pseudo Tim e Series p+d 0 2
Modeled Pseudo Tim e Series
q 1 0
Table 9. The models suggested by the SCAN option from PROC ARIMA. Models suggested by the ARIMA Procedure with ARMA (0,1) chosen as the applicable format. MA(1): (Yt t 1 t 1 ) is stationary for all for all ) (Hamilton, 1994). The following data set shows the Euclidean Distance values for the two curves and the index point . Index
Euclidean Distance for Model Curve
Euclidean Distance for Experim ental Curve
1 2 3 4 5 6 7 8 9 10 11
123.991 114.005 103.132 75.792 53.333 40.632 49.020 73.378 99.871 126.127 148.660
112.603 103.016 90.733 59.081 53.794 37.162 42.343 66.215 69.759 86.763 97.277
12 13 14 15 16 17 18 19 20 21 13
Euclidean Distance for Euclidean Distance for Model Curve Experim ental Curve 165.825 176.307 178.105 170.231 150.809 120.055 75.325 34.299 83.349 177.819 176.307
129.985 138.553 188.457 178.690 139.015 91.083 71.192 46.317 69.807 167.173 138.553
Table 10. Dataset Sum m ary of the Euclidean Distances for both curves. These w ere com pared using PROC Sim ilarity
The two pseudo time series are shown in an overlay in the graph below. It appears that they are in good agreement.
Figure 7. The tw o pseudo tim e series plotted on the sam e graph. The tim es series in red(input) is the experimental curve Pseudo Tim e Series. The tim e series in blue (target) is the m odeled Pseudo Tim e Series.
The table and warp map that follows show that the distances are kept relatively minimum. As can be seen from the dynamic time warp map, most of the deviation from similarity takes place from points 7 to 15. This is also noted in the original comparison of the experimental curve and the curve generated using control points for the cubic Bezier spline.
Figure 8. Path diagram for the minal paths for linking the nodes for both Pseudo Time Series.
Path Missing Map Direct Maps Com pression Expansion Warps
Num b. 0 18 3 3 6
Path % 0 75 12.5 12.5 25
Input % 0 85.71 14.29 14.29 28.57
Target % 0 85.71 14.29 14.29 28.57
Max 0 6 2 1 2
Path Maxim um Input Maxim um % % 0 0 25.00 28.57 8.333 9.524 4.167 4.762 8.333 9.524
Target Maxim um % 0 28.57 9.524 4.762 9.524
Table 11. This table show s that there w ere 6 w arps (3 com pressions and 3 expansion).
The warp diagram displays good agreement between the respective time series.
Figure 14. The trace in red is the experim ental curve and the trace in blue is the m odeled curve. This m ap show s the m inimum distance m ap for connecting nearest nodes. Note that the m ajority of the difference in distances com es from indices 5 to 15.
The minimum relative measure of connecting the nearest nodes of the respective time series is approximately 8.33. Minimum Measure Summary Input Variable Model Pseudo Time Series
Experimental Pseudo Time Series
Table 12. The m inim um m easure of connecting sequential elements of the tim e series.
Conclusion The processes involved in retro-fitting a Bezier-type curve to the experimentally plotted PK/PD curve produced a very similar model. The experimental curve in the present study is fairly easily modeled by a cubic Bezier spline. Non-simple curves and more complex experimental curves might require splines of higher order. PROC TRANSREG is very effective modeling splines, but, as found in this case, parametric equations are needed to model the curve as a relation and not as a function. The intent of Dynamic Time Warping is to find a minimally weighted sum of path lengths between the two compared series. In the present example, the sum of the Dynamic Time Warp distances is 8.32, a much smaller number than the distances plotted for the vertical axis of each of the two series. It is well understood that Dynamic Time Warping does not produce a distance metric in the usual sense in that it does not obey the triangle inequality. However, two series that are identical in magnitude and phase should, theoretically, produce a weighted Dynamic Time Warping of 0. Thus, a relatively small value compared to the magnitude of the series numerical values, would indicate a stronger similarity for the two series that are compared. A production of a model PK/PD curve that is very similar to the experimentally derived curve provides a chance to form a set of parametric equations that define the system under study. These parametrically defined equations can be composed with a time-parameter function to allow the prediction of affect response versus time. One could then calculate the percent response to a drug based upon the time elapsed since initial treatment. It is interesting that the solutions to the a priori modeling process (briefly described in Appendix 2) would produce equations similar to the composition functions of the Bezier equations with exponential time curves.
Appendices Appendix A Program adapted from Wicklin, R., Compute the centroid of a polygon in SAS, The DO Loop.
Appendix B Use of the Bouc-Wen model for hysteresis loop systems and the Hill function from Biochemistry are becoming more common in dose response curves, gene regulation networks, and studies in anesthesiology. Drugs and enzymes often synergistically interact in biological systems. Therefore, any drugs, their metabolites and activated enzymes must be accounted for in a system of equations and/or differential equations. The Hill function is: Equation 14
ax 0n b n x 0n
Where x0 may denote the drug concentration in the experiment. x 1 represents molecular structural changes in a hybrid form that allows the expression of subsequently expressed enzymes. The values for a and b are constants relevant to the system, while n is the Hill coefficient. The overall model for the system is designed from an adapted Bouc -Wen non-linear differential equation and its sub-functions.
dw dx dx dw dw and w x w dt dt dt dt dt Equation 16 n x1 x0 ax0min x 1 2 , c0 n n cm c0 b x0min
n ax0max x x , w 1 2 2 2min n n b x0max x2max x2min
Here x2 is the expressed enzyme concentration, w is the normalization of x2 and ρ, σ, γ, κ are model parameters determined through experiment.
Note: The above equations can be modeled by PROC MODEL.
Appendix C The following program shows the code for the Bezier Curve in Figure 5. Code suggested by Dr. Rick Wicklin (personal communication). data adapt; do tau = 0 to 1 by 0.05; x = 50 + 174*tau + 1689*tau**2 - 1913*tau**3; y = 58 + 219*tau - 555*tau**2 + 278*tau**3; output; end; run; proc sgplot data=adapt; series x=x y=y /markers datalabel=tau datalabelattrs=(size=16); run;
References Belal, S., Shashank, K. T., Shashank K., Md. Nadeem A., Satish, K.D., Rahul S. (2013) Inverse Point Solution of Bezier Curve, International Journal of Scientific & Engineering Research, 4(6) DeMedeiros-Martins, A., Doria-Neto, A.D., and DeNelo, J.D. (2004) Comparison between Mahalanobis distance and Kullback-Leibler divergence in clustering analysis. Available at http://www.wseas.us/elibrary/conferences/udine2004/papers/483-202.pdf Hamilton, J.D. (1994), Time Series Analysis, Princeton, NJ: Princeton University Press. Hu J., Qin, K. R., Xiang, C., Lee T. H. (2012), Modeling of Hysteresis in Gene Regulatory Networks, Bulletin of Mathematical Biology. 74(8) 1727-1753 Leonard, M., Sloan, J. Lee, T. and Alzheimer, B. (2008) “An Introduction to Similarity Analysis Using SAS”, Proceedings of International Symposium of Forecasting. SAS Institute Inc., Cary, NC. Available at: https://support.sas.com/rnd/app/ets/papers/similarityanalysis.pdf Louizos, C., Yanez, J., Forrest, L., Davies, N. (2014) Understanding the hysteresis Loop conundrum in pharmacokinetic / pharmacodynamic relationships. J. Pharm. Sci. 17(1) 34-91. Nygaard, N., Time Series (2006) Found at: www.finmath.uchicago.edu/docs/Time%20S eries.pdf Pleuvry. Barabara J. (2005) Hysteresis in Drug Response, Anaethesia & Intensive Care Medicine, 6(8), 286-287 Ramsay, J.O. and Silverman, B.W. (2005). Functional data analysis, 2 nd Ed. New York: Springer Rumpler, M.J., Kandala, B., Vickroy, T. W., Hochhaus, G., Samms, R. A., (2013) Pharmacokinetics and pharmacodynamics of glycopyrrolate following continuous-infusion in the horse. J Vet. Pharmacol. Therap. 37, 133-144 Taghhia, J, (2011) Forcing Bezier Interpolation. Available at: http://www.javadtaghia.com/control-systemmath/forcing-bezier-interpolation Tran, V., (2007) Reverse Engineering Bezier Curves. Found at: Polymath Programmer polymathprogrammer.com/2007/06/27/reverse-engineering-bezier-curves/ Wei L., E Keogh, X Xi, M Yoder, (2008) Efficiently finding unusual shapes in large image databases. Data Mining and Knowledge Discovery, New York: Springer Wicklin, R., (2012) Testing data for multivariate normality, SAS Blogs Home > The DO Loop Available at: http://blogs.sas.com/content/iml/2012/03/02/testing-data-for-multivariate-normality.html Wicklin, R., (2012) Compute the multivariate normal density in SAS, SAS Blogs Home > The DO Loop Available at: http://blogs.sas.com/content/iml/2012/07/05/compute-multivariate-normal-denstity.html Wicklin, R., (2016) Compute the centroid of a polygon in SAS, SAS Blogs Home > The DO Loop Available at: http://blogs.sas.com/content/iml/2016/01/13/compute-centroid-polygon-sas.html Zhang, D. and Lu, G. (2004). Review of shape representation and description techniques. Pattern Recognition, 37(1): 1-19.
Acknowledgements The author would like to thank Rick Wicklin, PhD, Distinguished Researcher at the SAS Institute in Cary, North Carolina and Arthur Li, Professor and Researcher at the City of Hope, Duarte, California. These gentlemen provided much needed and much appreciated structural and editorial support for this paper.
Contact Information Your comments and questions are valued and encouraged. Contact the author at: Ronald Smith Florida Southwestern State College / Professor and Statistical Consultant (239) 292-1810 [email protected]
SAS and all other SAS Institute Inc. product or service names are registered trademarks or trademarks of SAS Institute Inc. in the USA and other countries. ®indicates USA registration. Other brand and product names are trademarks of their respective companies.