- Email: [email protected]
Mathematical model and computer representation for nonlinear response systems
COMPUTERS
AND
BIOMEDICAL
Mathematical
RESEARCH
Model Nonlinear
(1972)
5,239-246
and Computer Representation Response Systems
for
W. B. QUAY, N. G. SMIRIGA, AND W. S. HAUGELAND Department of Zoology, University of California, Berkeley, California 94720 Department of Mathematics, San Francisco State College, San Francisco, California 94132 Department of Computer Science, University of California, Berkeley, California 94720
ReceivedOctober 6.1971 Thebeta-distribution isproposedasa suitablyflexiblecumulativedistributionfunction for the representation of data andthe facilitation of analysisfor a widevariety of nonlinearresponse systemsshowingadaptation,or approachingsaturation.A particularly difficult response systeminvolving behavioralphase-shifting followingreversalin timing of photoperiod(Zeitgeber)is usedasprimary test subject.The modelfor this system becomes Ix,&,
p) =
j&l2
- uy-’ du
II
12”+J-‘B(a,
/3),
where0 G X G 12.Computational variablesandlimitationsarediscussed andappropriate computersubprograms arenoted.
A common and difficult problem in studies on complex regulatory mechanisms is the existing deficiency in adequate and flexible mathematical models for the representation of nonlinear responses. Such models are needed for optimal representation by curves, of empirical and variable data from nonlinear response systems. Of particular concern are those classesof responsesthat are either wholly or partially non-
linear, whatever category of graphic scale or transformation is applied to them. A case of this sort prompted the study upon which this report is based. This particular case concerns neural and neuroendoerine regulation of experimental phase-shifts in mammalian 24-hour rhythms. This can be considered to be a complex biological stimulus-response systemshowing adaptation. Sensitivity in such a system changes with stimulus intensity, but the change occurs after a lag, and kinetics of adaptation must be analyzed in order to understand stimulus-responseevents (I). Responseor adaptation curves obtained in this work resemble those from a wide variety of phenomena, such as water absorption by proteins in different water vapor pressures 0 1972 by Academic
Press,
Inc.
239
240
QUAY, SMIRIGA, AND HAUGELAND
(2), mortality after different doses of radiation (3), enzyme activity in the presence of allosteric activators (4), and synaptic input-output relations in neurons (5). The primary aim of this study was the development of a broadly flexible model for stimulus-response and adaptation phenomena, with particular emphasis on possibilities for the derivation from it of rigorous analytical and statistical methods. The model and the derivative methods are considered to be most useful in the following ways: (1) to allow construction of a definitive (unique answer) response or adaptation curve for variable experimental data from either single individuals or series of pooled ones; (2) to enable exact comparisons to be made between two or more curves derived in this way, and the analysis of their variances and differences by statistical methods; (3) to establish a macroscopic basis or scale upon which the revealing and analysis of microscopic phasic features will be possible with suitable mathematical reliability and sensitivity; (4) to facilitate the development of theories for either particular or categorical nonlinear response systems; and (5) to lead to an efficient and quantitative modus operandi for experimentally testing such theories. The use of this model and the derivative methods is excessively laborious without employing a computer. In this study, development and testing of the mathematical model was accompanied by the development and modification of suitable computer subprograms for use in the CDC 6400. EXAMPLENONLINEARRESPONSESYSTEM
The nonlinear response system used as a test subject and example in this study is the phase-shifting of the circadian or 24-hour rhythm in running activity by laboratory rats. In man and other mammals there are more than 100 functions and structural elements which have such a daily rhythm (6), and which phase-shift when the timing of environmental variables or the daily schedule is changed. The rhythm used here, that of spontaneous running activity by female rats, can be considered representative. The most sharply timed characteristic in this rhythm in our experimental conditions is usually the daily start of major running activity. Onset of darkness is most closely correlated with this event and may be considered to be the cue or Zeitgeber for the response and timing of running activity. The daily period of light or photoperiod, fixed in duration and intensity, is either the dominant Zeitgeber or is directly related to it. Experimental daily light and dark periods are 12 hours in length (LD 12: 12) and are precisely timed automatically. The daily schedule of the photoperiod is ordinarily fixed, but when it is shifted (AD) or reversed (A@ = 12) the start of running activity shifts also (d&S). 04s is a nonlinear and individually variable response, requiring about a week for completion. Completion of d&S is entrainment. Graphs and data points from an example experiment (7) are shown in Fig. 1. THE MODEL
Y =f(X) denote the relationship between A’ (-days post A@) and Y (=A@= % sh’ft I in starting time). Then, Y =,f(X) is a nondecreasing function for Let
NONLINEAR
RESPONSE
241
SYSTEMS
which 0 G Y G 1. Therefore, a cumulative distribution function should be a good model for the experimental data for O&S (Fig. 1). Fundamentally, a cumulative distribution function is nonnegative, nondecreasing and F(-co) = 0, F(+a) = 1. Moreover, this cumulative distribution should be very flexible so that observed extremes of &B(Fig. 1) would give a good fit. Such extremes may be typified by those
0
I
2
3
4
5
6
7
*+SP,
observed computed AL-- __ __ . _. . . . .
&Sip
.
A+Sxs+,
. -----
*$&,
o -
8
9
-----
IO
II
I2
X (= days post A@, FIG. 1. Results of an experiment from the nonlinear responsesystemused as an example in this study. Observed data points and computed curves are shown for extreme individual pinealectomized (d&S’,,) and sham-operated (d#?,,) animals and for means of pooled data for pinealectomized (A&,) and control (d&S,,+,) aniamls. Further explanation in text. Data from (7).
individual pinealectomized (Fig. 1, A$$,,) and sham-operated (Fig. 1, A$&,) animals whose d&S’ data and curves are figured here. The beta-distribution offers the flexibility required for the observed data. Recall that the beta-distribution is, by definition, given by Lx% PI = &(a, B)/~(a~ B) (1) where B,(a,p) is the incomplete
beta-function,
i.e.,
&A% B> = j: ZP-‘(1 - .?+-I d’u
(2)
0
and B(cr,/3) is the beta-function,
i.e.,
Nay B> = ~<4WW~
+ 6%
(3)
242
QUAY,SMIRIGA,
where r is the gamma-function,
AND HAUGELAND
i.e.,
r(a) = Jua-1 e-udu.
(4)
0
Further information on the beta-distribution is available (8). In this definition the domain of X is 0 G X < I. Since the domain of the example experimental data is 0 G XG 12, where A&? is studied over 12 days after d@, the model for these data becomes
where 0 f Xd 12. COMPUTATIONS
A nonlinear regression program is used to fit the data to the distribution I,,,,(oz,/~). This means that the program attempts to find those values of LYand fl which will minimize the relative sum of squares
where Xj =jis thej-th day, and Yj is the percentage ofthe shift in starting time (d&S) observed at time Xj. As mentioned before, the beta-distribution has been chosen because of its flexibility in shape; however, the evaluation of Ixjl12(cc, ,!?)is somewhat tedious. The fit has been done by a modification of a previously described (9) nonlinear regression program. For this published program the user has to supply zxj/12(c(,
B> z Bxjl12(c(3P)/B(a, P).
(7)
The integral Bxj, iZ(c(,p) has been numerically evaluated using Wynn’s (10) algorithm and the beta-function B(cr,/3) = r(cc)r@)/I’(a + /3) has been evaluated by using Wrench’s (10) algorithm. The numerical evaluation of Ix,, ,z(cc,p) thus obtained has been compared for several values of X,i, ccand /3 with that obtained through use of Pearson’s (8) tables. These tables provide only seven digits after the decimal point. The computations yield thirteen digits (when done on a CDC 6400). Rounding the results to seven digits showed a perfect agreement with the tables. The least-squares program requires as input an initial guess for the parameters ccand p. The Biometrika Tables (12) are very convenient for this. Actually, however, the choice of the initial values of CLand p is not too important because the program is not very sensitive to the initial choice of the parameters; convergence is achieved rapidly and in most cases occurs before twenty iterations (Table 1).
NONLINEAR
RESPONSE
243
SYSTEMS
Characteristics of the input data for X and Y are important in the design of computations and analyses that will derive the best fit and that will be most instructive concerning the composition of complex nonlinear response systems. In the case of the response system used as a primary example here (Fig. l), Xis an independent and TABLE
1
EFFECTS OF INITIAL GUESSES FOR THE PARAMETERS (Y AND /3 ON COMPUTED VALUES, USING DATA FROM THE EXPERIMENT ILLUSTRATED BY FIG. 1. “Guesses” Sample
Computed B
2.750 5.200 2.000 7.400 5.900
2.590 8.300 2.500 7.900 6.880
1.142436 1.142439 1.142439 1.142438 1.142434
1.388895 1.388895 1.388897 1.388898 1.388896
i
.090817
t
.090817 .090816
$
.090816
1.500 2.000 2.500 5.001 4.010
3.500 4.500 6.000 12.001 8.010
1.027620 1.027615 1.027616 1.027617 1.027615
2.240555 2.240543 2.240546 2.240548 2.240542
:
.016202 .016202
i 4
.016202 .016202
Ah%,,
1.500 1.190 1.500 .900 2.500
15.000 3.580 9.600 3.200 30.000
1.039214 .2040684 1.022777 .2040685 2.110454
12.12654 1.813353 12.37003 1.813354 23.66457
It 4 L 4 L
.053175 .010262 .053908 .010262 1.9553
4%
1.100 1.000 .500 .900 .750
15.000 4.500 3.500 12.560 3.000
1.074941 .001119 .004167066 .8478810 .004167079
13.46208 .3760267 .1225355 13.81600 .1225359
It L 4 i 4
AWs,
114 = reduced; % = termination L = limited by preset iterations
a
Sum of squares effect” final sum
a
due to neglible (=40).
P
change
in sum
of squares
.10219 .011420 SW12857 .052135 .0012857 (nonconvergence);
fixed variable. Standardization of the number of X values employed in an experimental design is advantageous for efficient use of the computer program and the comparisons of replications of diverse kinds. In the case of studies on phase shifts, such as those in Fig. 1, the number of values of X (days) should exceed the number necessary for unequivocal saturation (entrainment) to occur. Although the number of values of X necessary for saturation may be easily apparent in many studies; in others, such as one recently reported (13), this is not always true. The degree to which saturation is real or approximated thus can vary.
244
QUAY,SMIRIGA,
AND HAUGELAND
The fact that Xis an independent and fixed variable while Y is dependent and continuous in the case of the primary example system used here (Fig. 1) leads to certain limitations in options for statistical analyses of the data. But, this of itself is not limiting in relation to the applicability of the model to cumulative distribution functions in which the variables X and Y are both continuous or may have certain other properties. The continuous variable Y can have scale and reference points that are either directly observational or that are some function of observed values, as in the example of Fig. 1. Here the time span (12 hours) for complete shift in activity starting time (d@) follows from the change (12 hours) in the new time of the onset of darkness after the reversal or shift in daily timing of the environmental photoperiod. In some experimental situations, however, the entrainment endpoint, or point of saturation for fl@, may become dissociated from the Zeitgeber endpoint (O$S # 443) (13). In this case, computational application of the model requires the calculation of a new and arbitrary transformation of the observational Y scale in hours to values of O-l .O. In our primary example situation (Fig. l), nevertheless, Y = 1.O. or saturation, when O$S = 12 hours = d@, =lOO% shift in the timing of the response. Assignment of values of Y to particular data does not depend only on selection or transformation of the Y scale as described above. Assignment of Y values depends, also, on the exact determination of the points of origin ( Y = 0) and saturation ( Y ~1.O) for a particular set of data. This requires having either invariable and predictable points of origin and saturation, or sufficient replicate measurements at these points to give them validity. It so happens that in the response system illustrated in Fig. I. the observational points for Y = 0 and Y = 1.O in O$S are often slightly but significantly before or after those for d@. Observational data points for Y may not be of uniform quality or validity, at least in terms of their being uniform in statistical variance. Many experiments with the response system in Fig. 1 have shown that the coefficiznt of variation for the sequence of observational Y data points (from X = 1 to X = 12) follows a nonlinear curve. Coefficients of variance or other measures of variability are calculated from pooled data points or from replications with mutliple test animals or with multiple phase shifts. The data points for Y, as input to the computer program, are then weighted inversely in relation to their variance. RESULTS AND DISCUSSION
Our model has accommodated data plots from a variety of experimental manipulations of the primary response system under study. One set of these plots of observed and computed points is shown in Fig. I. Another set is illustrated here by the computed curves only (Fig. 2). In this experiment a transitory difference between O#JS,; and Or&; is illustrated and is substantiated by statistical treatment of the
245
NONLINEAR RESPONSE SYSTEMS
@ 3.2 days
0
I
2
3
4
5
6
7
8
9
10 II
I;
X
I
2
3
4
5
6 X
7
8
9
IO
II
r I2
FIG. 2. Experimental results from the nonlinear response system used as an example in this study. Computed response curves of matched samples of siblings with the pineal stalk cut (O#S,J or with a sham-operation (d+Q are based on means inversely weighted with the coefficients of variation. Further explanation in text. Data from (13).
data (13). These and other related results with this system suggest that effective comparative analyses can be made of separate segments of the response curves, in attempts to distinguish the temporal relations and relative importance of different sources of regulatory input. The potential resolving power of such attempts is yet to be fully explored but can be expected to vary with the particular segment of the continuous variable Y. APPENDIX
Outline of subroutines modified and developed for the computer (FORTRAN, CDC 6400program and used in this study. Numbers l-5 were modified from Baer (9). 1. Subroutine NLIN: PROGRAM NLINEAR (INPUT, OUTPUT) 2. Subroutine INDATA : SUBROUTINE INDATA (CNTROL, B, Z, MDIM, NDIM, BDIM) 3. Subroutine SETUP : SUBROUTINE SETUP (CNTROL, B, Z, A, SVMSQ, MDIM, BDIM) 4. Subroutine GAUSS SUBROUTINE GAUSS (CNTROL, B. Z, MDIM, NDIM, BDIM, A) 5. Subroutine FINALE: SUBROUTINE FINALE (CNTROL, B, Z, MDIM, NDIM) 6. Subroutine DMPMAT: SUBROUTINE DMPMAT (CNTROL, M,A,C)
246
QUAY, SMIRIGA, AND HAUGELAND
7. Subroutine DERIV: FUNCTION DERIV (K, N, B, Z, MDIM) 8. Subroutine YCOMP: FUNCTION YCOMP (N, B, Z) 9. Subroutine CCF: FUNCTION CCF (X, P, Q) 10. Subroutine GAMM FUNCTION GAMM (X) ACKNOWLEDGMENTS We are grateful to Professors L. LeCam (Department of Statistics, University of California, Berkeley) and B. N. Parlett (Department of Computer Science, University of California, Berkeley) for important discussions and advice during this study. This work was supported by a research grant (NS-06296) from the National institutes of Health, U. S. Public Health Service, and aid from Biomedical Sciences Support Grant funds, and the Miller Institute for Basic Research in Science, University of California, Berkeley.
REFERENCES 1. VARJ$ D. Uber Nichtlineare Analogschaltungen zur Simulierung Biologischer Adaptationsvorglnge. Progr. Bruin Res. 17, 74 (1965). 2. STEINHARDT, J. AND REYNOLDS,J. A. “MultipleEquilibriainProteins,“p. 154. Academic Press, New York, 1969. 3. IBERALL, A. S. Quantitative modeling of the physiological factors in radiation lethality. Ann. N. Y. Acad. Sci. 147, 1 (1967). 4. HESS, B. Biochemical regulations. In “Systems Theory and Biology: Proceedings of the III Systems Symposium at Case Institute of Technology” (M. D. Mesarovic, Ed.), p. 88. SpringerVerlag, New York, 1968. 5. TAYLOR, W. K. A model of learning mechanisms in the brain. Progr. Bruin Res. 17,369 (1965). 6. ASCHOFF,J. Desynchronization and resynchronization of human circadian rhythms. Aerospnce Med. 60, 844 (1969). 7. QUAY, W. B. Precocious entrainment and associated characteristics of activity patterns following pinealectomy and reversal of photoperiod. Physiu[. Beh~v. 5, 1281 (1970). 8. PEARSON, K. “Tables of the Incomplete Beta-function,” 2nd ed., Cambridge Univ. Press, Cambridge, 1968. 9. BAER, R. M. Nonlinear regression and the solution of simultaneous equations. Comm. Assoc. Comput. Mach. 5, 397 (1962). 10. WYNN, P. “An Arsenal of Algol Procedures for the Evaluation of Continuous Fractions and for Effecting the Epsilon Algorithm,” p. 36. MRC Technical Summary Report No. 537, Mathematics Research Center, United States Army, University of Wisconsin, Madison, WI, 1965. Il. WRENCH, J. W., Jr. Concerning two series for the gamma function. Math. Comput. 22, 617 (1968). 12. PEARSON,E. S. AND HARTLEY, H. 0. “Biometrika Tables for Statisticians,” Cambridge Univ. Press, Cambridge, 1954. 13. QUAY, W. B. Dissimilar functional effects of pineal stalk and cerebral meningeal interruptions on phase shifts of circadian activity rhythms. Physiul. Be&v. 7, 557 (1971).
AND
BIOMEDICAL
Mathematical
RESEARCH
Model Nonlinear
(1972)
5,239-246
and Computer Representation Response Systems
for
W. B. QUAY, N. G. SMIRIGA, AND W. S. HAUGELAND Department of Zoology, University of California, Berkeley, California 94720 Department of Mathematics, San Francisco State College, San Francisco, California 94132 Department of Computer Science, University of California, Berkeley, California 94720
ReceivedOctober 6.1971 Thebeta-distribution isproposedasa suitablyflexiblecumulativedistributionfunction for the representation of data andthe facilitation of analysisfor a widevariety of nonlinearresponse systemsshowingadaptation,or approachingsaturation.A particularly difficult response systeminvolving behavioralphase-shifting followingreversalin timing of photoperiod(Zeitgeber)is usedasprimary test subject.The modelfor this system becomes Ix,&,
p) =
j&l2
- uy-’ du
II
12”+J-‘B(a,
/3),
where0 G X G 12.Computational variablesandlimitationsarediscussed andappropriate computersubprograms arenoted.
A common and difficult problem in studies on complex regulatory mechanisms is the existing deficiency in adequate and flexible mathematical models for the representation of nonlinear responses. Such models are needed for optimal representation by curves, of empirical and variable data from nonlinear response systems. Of particular concern are those classesof responsesthat are either wholly or partially non-
linear, whatever category of graphic scale or transformation is applied to them. A case of this sort prompted the study upon which this report is based. This particular case concerns neural and neuroendoerine regulation of experimental phase-shifts in mammalian 24-hour rhythms. This can be considered to be a complex biological stimulus-response systemshowing adaptation. Sensitivity in such a system changes with stimulus intensity, but the change occurs after a lag, and kinetics of adaptation must be analyzed in order to understand stimulus-responseevents (I). Responseor adaptation curves obtained in this work resemble those from a wide variety of phenomena, such as water absorption by proteins in different water vapor pressures 0 1972 by Academic
Press,
Inc.
239
240
QUAY, SMIRIGA, AND HAUGELAND
(2), mortality after different doses of radiation (3), enzyme activity in the presence of allosteric activators (4), and synaptic input-output relations in neurons (5). The primary aim of this study was the development of a broadly flexible model for stimulus-response and adaptation phenomena, with particular emphasis on possibilities for the derivation from it of rigorous analytical and statistical methods. The model and the derivative methods are considered to be most useful in the following ways: (1) to allow construction of a definitive (unique answer) response or adaptation curve for variable experimental data from either single individuals or series of pooled ones; (2) to enable exact comparisons to be made between two or more curves derived in this way, and the analysis of their variances and differences by statistical methods; (3) to establish a macroscopic basis or scale upon which the revealing and analysis of microscopic phasic features will be possible with suitable mathematical reliability and sensitivity; (4) to facilitate the development of theories for either particular or categorical nonlinear response systems; and (5) to lead to an efficient and quantitative modus operandi for experimentally testing such theories. The use of this model and the derivative methods is excessively laborious without employing a computer. In this study, development and testing of the mathematical model was accompanied by the development and modification of suitable computer subprograms for use in the CDC 6400. EXAMPLENONLINEARRESPONSESYSTEM
The nonlinear response system used as a test subject and example in this study is the phase-shifting of the circadian or 24-hour rhythm in running activity by laboratory rats. In man and other mammals there are more than 100 functions and structural elements which have such a daily rhythm (6), and which phase-shift when the timing of environmental variables or the daily schedule is changed. The rhythm used here, that of spontaneous running activity by female rats, can be considered representative. The most sharply timed characteristic in this rhythm in our experimental conditions is usually the daily start of major running activity. Onset of darkness is most closely correlated with this event and may be considered to be the cue or Zeitgeber for the response and timing of running activity. The daily period of light or photoperiod, fixed in duration and intensity, is either the dominant Zeitgeber or is directly related to it. Experimental daily light and dark periods are 12 hours in length (LD 12: 12) and are precisely timed automatically. The daily schedule of the photoperiod is ordinarily fixed, but when it is shifted (AD) or reversed (A@ = 12) the start of running activity shifts also (d&S). 04s is a nonlinear and individually variable response, requiring about a week for completion. Completion of d&S is entrainment. Graphs and data points from an example experiment (7) are shown in Fig. 1. THE MODEL
Y =f(X) denote the relationship between A’ (-days post A@) and Y (=A@= % sh’ft I in starting time). Then, Y =,f(X) is a nondecreasing function for Let
NONLINEAR
RESPONSE
241
SYSTEMS
which 0 G Y G 1. Therefore, a cumulative distribution function should be a good model for the experimental data for O&S (Fig. 1). Fundamentally, a cumulative distribution function is nonnegative, nondecreasing and F(-co) = 0, F(+a) = 1. Moreover, this cumulative distribution should be very flexible so that observed extremes of &B(Fig. 1) would give a good fit. Such extremes may be typified by those
0
I
2
3
4
5
6
7
*+SP,
observed computed AL-- __ __ . _. . . . .
&Sip
.
A+Sxs+,
. -----
*$&,
o -
8
9
-----
IO
II
I2
X (= days post A@, FIG. 1. Results of an experiment from the nonlinear responsesystemused as an example in this study. Observed data points and computed curves are shown for extreme individual pinealectomized (d&S’,,) and sham-operated (d#?,,) animals and for means of pooled data for pinealectomized (A&,) and control (d&S,,+,) aniamls. Further explanation in text. Data from (7).
individual pinealectomized (Fig. 1, A$$,,) and sham-operated (Fig. 1, A$&,) animals whose d&S’ data and curves are figured here. The beta-distribution offers the flexibility required for the observed data. Recall that the beta-distribution is, by definition, given by Lx% PI = &(a, B)/~(a~ B) (1) where B,(a,p) is the incomplete
beta-function,
i.e.,
&A% B> = j: ZP-‘(1 - .?+-I d’u
(2)
0
and B(cr,/3) is the beta-function,
i.e.,
Nay B> = ~<4WW~
+ 6%
(3)
242
QUAY,SMIRIGA,
where r is the gamma-function,
AND HAUGELAND
i.e.,
r(a) = Jua-1 e-udu.
(4)
0
Further information on the beta-distribution is available (8). In this definition the domain of X is 0 G X < I. Since the domain of the example experimental data is 0 G XG 12, where A&? is studied over 12 days after d@, the model for these data becomes
where 0 f Xd 12. COMPUTATIONS
A nonlinear regression program is used to fit the data to the distribution I,,,,(oz,/~). This means that the program attempts to find those values of LYand fl which will minimize the relative sum of squares
where Xj =jis thej-th day, and Yj is the percentage ofthe shift in starting time (d&S) observed at time Xj. As mentioned before, the beta-distribution has been chosen because of its flexibility in shape; however, the evaluation of Ixjl12(cc, ,!?)is somewhat tedious. The fit has been done by a modification of a previously described (9) nonlinear regression program. For this published program the user has to supply zxj/12(c(,
B> z Bxjl12(c(3P)/B(a, P).
(7)
The integral Bxj, iZ(c(,p) has been numerically evaluated using Wynn’s (10) algorithm and the beta-function B(cr,/3) = r(cc)r@)/I’(a + /3) has been evaluated by using Wrench’s (10) algorithm. The numerical evaluation of Ix,, ,z(cc,p) thus obtained has been compared for several values of X,i, ccand /3 with that obtained through use of Pearson’s (8) tables. These tables provide only seven digits after the decimal point. The computations yield thirteen digits (when done on a CDC 6400). Rounding the results to seven digits showed a perfect agreement with the tables. The least-squares program requires as input an initial guess for the parameters ccand p. The Biometrika Tables (12) are very convenient for this. Actually, however, the choice of the initial values of CLand p is not too important because the program is not very sensitive to the initial choice of the parameters; convergence is achieved rapidly and in most cases occurs before twenty iterations (Table 1).
NONLINEAR
RESPONSE
243
SYSTEMS
Characteristics of the input data for X and Y are important in the design of computations and analyses that will derive the best fit and that will be most instructive concerning the composition of complex nonlinear response systems. In the case of the response system used as a primary example here (Fig. l), Xis an independent and TABLE
1
EFFECTS OF INITIAL GUESSES FOR THE PARAMETERS (Y AND /3 ON COMPUTED VALUES, USING DATA FROM THE EXPERIMENT ILLUSTRATED BY FIG. 1. “Guesses” Sample
Computed B
2.750 5.200 2.000 7.400 5.900
2.590 8.300 2.500 7.900 6.880
1.142436 1.142439 1.142439 1.142438 1.142434
1.388895 1.388895 1.388897 1.388898 1.388896
i
.090817
t
.090817 .090816
$
.090816
1.500 2.000 2.500 5.001 4.010
3.500 4.500 6.000 12.001 8.010
1.027620 1.027615 1.027616 1.027617 1.027615
2.240555 2.240543 2.240546 2.240548 2.240542
:
.016202 .016202
i 4
.016202 .016202
Ah%,,
1.500 1.190 1.500 .900 2.500
15.000 3.580 9.600 3.200 30.000
1.039214 .2040684 1.022777 .2040685 2.110454
12.12654 1.813353 12.37003 1.813354 23.66457
It 4 L 4 L
.053175 .010262 .053908 .010262 1.9553
4%
1.100 1.000 .500 .900 .750
15.000 4.500 3.500 12.560 3.000
1.074941 .001119 .004167066 .8478810 .004167079
13.46208 .3760267 .1225355 13.81600 .1225359
It L 4 i 4
AWs,
114 = reduced; % = termination L = limited by preset iterations
a
Sum of squares effect” final sum
a
due to neglible (=40).
P
change
in sum
of squares
.10219 .011420 SW12857 .052135 .0012857 (nonconvergence);
fixed variable. Standardization of the number of X values employed in an experimental design is advantageous for efficient use of the computer program and the comparisons of replications of diverse kinds. In the case of studies on phase shifts, such as those in Fig. 1, the number of values of X (days) should exceed the number necessary for unequivocal saturation (entrainment) to occur. Although the number of values of X necessary for saturation may be easily apparent in many studies; in others, such as one recently reported (13), this is not always true. The degree to which saturation is real or approximated thus can vary.
244
QUAY,SMIRIGA,
AND HAUGELAND
The fact that Xis an independent and fixed variable while Y is dependent and continuous in the case of the primary example system used here (Fig. 1) leads to certain limitations in options for statistical analyses of the data. But, this of itself is not limiting in relation to the applicability of the model to cumulative distribution functions in which the variables X and Y are both continuous or may have certain other properties. The continuous variable Y can have scale and reference points that are either directly observational or that are some function of observed values, as in the example of Fig. 1. Here the time span (12 hours) for complete shift in activity starting time (d@) follows from the change (12 hours) in the new time of the onset of darkness after the reversal or shift in daily timing of the environmental photoperiod. In some experimental situations, however, the entrainment endpoint, or point of saturation for fl@, may become dissociated from the Zeitgeber endpoint (O$S # 443) (13). In this case, computational application of the model requires the calculation of a new and arbitrary transformation of the observational Y scale in hours to values of O-l .O. In our primary example situation (Fig. l), nevertheless, Y = 1.O. or saturation, when O$S = 12 hours = d@, =lOO% shift in the timing of the response. Assignment of values of Y to particular data does not depend only on selection or transformation of the Y scale as described above. Assignment of Y values depends, also, on the exact determination of the points of origin ( Y = 0) and saturation ( Y ~1.O) for a particular set of data. This requires having either invariable and predictable points of origin and saturation, or sufficient replicate measurements at these points to give them validity. It so happens that in the response system illustrated in Fig. I. the observational points for Y = 0 and Y = 1.O in O$S are often slightly but significantly before or after those for d@. Observational data points for Y may not be of uniform quality or validity, at least in terms of their being uniform in statistical variance. Many experiments with the response system in Fig. 1 have shown that the coefficiznt of variation for the sequence of observational Y data points (from X = 1 to X = 12) follows a nonlinear curve. Coefficients of variance or other measures of variability are calculated from pooled data points or from replications with mutliple test animals or with multiple phase shifts. The data points for Y, as input to the computer program, are then weighted inversely in relation to their variance. RESULTS AND DISCUSSION
Our model has accommodated data plots from a variety of experimental manipulations of the primary response system under study. One set of these plots of observed and computed points is shown in Fig. I. Another set is illustrated here by the computed curves only (Fig. 2). In this experiment a transitory difference between O#JS,; and Or&; is illustrated and is substantiated by statistical treatment of the
245
NONLINEAR RESPONSE SYSTEMS
@ 3.2 days
0
I
2
3
4
5
6
7
8
9
10 II
I;
X
I
2
3
4
5
6 X
7
8
9
IO
II
r I2
FIG. 2. Experimental results from the nonlinear response system used as an example in this study. Computed response curves of matched samples of siblings with the pineal stalk cut (O#S,J or with a sham-operation (d+Q are based on means inversely weighted with the coefficients of variation. Further explanation in text. Data from (13).
data (13). These and other related results with this system suggest that effective comparative analyses can be made of separate segments of the response curves, in attempts to distinguish the temporal relations and relative importance of different sources of regulatory input. The potential resolving power of such attempts is yet to be fully explored but can be expected to vary with the particular segment of the continuous variable Y. APPENDIX
Outline of subroutines modified and developed for the computer (FORTRAN, CDC 6400program and used in this study. Numbers l-5 were modified from Baer (9). 1. Subroutine NLIN: PROGRAM NLINEAR (INPUT, OUTPUT) 2. Subroutine INDATA : SUBROUTINE INDATA (CNTROL, B, Z, MDIM, NDIM, BDIM) 3. Subroutine SETUP : SUBROUTINE SETUP (CNTROL, B, Z, A, SVMSQ, MDIM, BDIM) 4. Subroutine GAUSS SUBROUTINE GAUSS (CNTROL, B. Z, MDIM, NDIM, BDIM, A) 5. Subroutine FINALE: SUBROUTINE FINALE (CNTROL, B, Z, MDIM, NDIM) 6. Subroutine DMPMAT: SUBROUTINE DMPMAT (CNTROL, M,A,C)
246
QUAY, SMIRIGA, AND HAUGELAND
7. Subroutine DERIV: FUNCTION DERIV (K, N, B, Z, MDIM) 8. Subroutine YCOMP: FUNCTION YCOMP (N, B, Z) 9. Subroutine CCF: FUNCTION CCF (X, P, Q) 10. Subroutine GAMM FUNCTION GAMM (X) ACKNOWLEDGMENTS We are grateful to Professors L. LeCam (Department of Statistics, University of California, Berkeley) and B. N. Parlett (Department of Computer Science, University of California, Berkeley) for important discussions and advice during this study. This work was supported by a research grant (NS-06296) from the National institutes of Health, U. S. Public Health Service, and aid from Biomedical Sciences Support Grant funds, and the Miller Institute for Basic Research in Science, University of California, Berkeley.
REFERENCES 1. VARJ$ D. Uber Nichtlineare Analogschaltungen zur Simulierung Biologischer Adaptationsvorglnge. Progr. Bruin Res. 17, 74 (1965). 2. STEINHARDT, J. AND REYNOLDS,J. A. “MultipleEquilibriainProteins,“p. 154. Academic Press, New York, 1969. 3. IBERALL, A. S. Quantitative modeling of the physiological factors in radiation lethality. Ann. N. Y. Acad. Sci. 147, 1 (1967). 4. HESS, B. Biochemical regulations. In “Systems Theory and Biology: Proceedings of the III Systems Symposium at Case Institute of Technology” (M. D. Mesarovic, Ed.), p. 88. SpringerVerlag, New York, 1968. 5. TAYLOR, W. K. A model of learning mechanisms in the brain. Progr. Bruin Res. 17,369 (1965). 6. ASCHOFF,J. Desynchronization and resynchronization of human circadian rhythms. Aerospnce Med. 60, 844 (1969). 7. QUAY, W. B. Precocious entrainment and associated characteristics of activity patterns following pinealectomy and reversal of photoperiod. Physiu[. Beh~v. 5, 1281 (1970). 8. PEARSON, K. “Tables of the Incomplete Beta-function,” 2nd ed., Cambridge Univ. Press, Cambridge, 1968. 9. BAER, R. M. Nonlinear regression and the solution of simultaneous equations. Comm. Assoc. Comput. Mach. 5, 397 (1962). 10. WYNN, P. “An Arsenal of Algol Procedures for the Evaluation of Continuous Fractions and for Effecting the Epsilon Algorithm,” p. 36. MRC Technical Summary Report No. 537, Mathematics Research Center, United States Army, University of Wisconsin, Madison, WI, 1965. Il. WRENCH, J. W., Jr. Concerning two series for the gamma function. Math. Comput. 22, 617 (1968). 12. PEARSON,E. S. AND HARTLEY, H. 0. “Biometrika Tables for Statisticians,” Cambridge Univ. Press, Cambridge, 1954. 13. QUAY, W. B. Dissimilar functional effects of pineal stalk and cerebral meningeal interruptions on phase shifts of circadian activity rhythms. Physiul. Be&v. 7, 557 (1971).