UCRL-TR-208164

Model-based Processing of
Micro-cantilever Sensor Arrays

J. W. Tringe, D. S. Clague, J. V. Candy, C. L. Lee,

R. E. Rudd, A. K. Burnham

November 23, 2004

Disclaimer

This document was prepared as an account of work sponsored by an agency of the United States
Government. Neither the United States Government nor the University of California nor any of their
employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for
the accuracy, completeness, or usefulness of any information, apparatus, product, or process
disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise,
does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United
States Government or the University of California. The views and opinions of authors expressed herein
do not necessarily state or reflect those of the United States Government or the University of California,
and shall not be used for advertising or product endorsement purposes.

This work was performed under the auspices of the U.S. Department of Energy by University of
California, Lawrence Livermore National Laboratory under Contract W-7405-Eng-48.
 MODEL-BASED PROCESSING OF MICRO-CANTILEVER SENSOR ARRAYS

J. Tringe, D. Clague, J. Candy, C. Lee, R. Rudd, and A. Burnham

University of California
Lawrence Livermore National Laboratory
Livermore, CA 94550

Executive Summary

We develop a model-based processor (MBP) for a micro-cantilever array sensor to detect target
species in solution. After discussing the generalized framework for this problem, we develop the specific
model used in this study. We perform a proof-of-concept experiment, fit the model parameters to the
measured data and use them to develop a Gauss-Markov simulation. We then investigate two cases of
interest: (1) averaged deflection data, and (2) multi-channel data. In both cases the evaluation proceeds by
first performing a model-based parameter estimation to extract the model parameters, next performing a
Gauss-Markov simulation, designing the optimal MBP and finally applying it to measured experimental
data. The simulation is used to evaluate the performance of the MBP in the multi-channel case and compare
it to a smoother (averager) typically used in this application. It was shown that the MBP not only
provides a significant gain (~ 80dB) in signal-to-noise ratio (SNR), but also consistently outperforms the
smoother by 40-60 dB. Finally, we apply the processor to the smoothed experimental data and demonstrate
its capability for chemical detection. The MBP performs quite well, though it includes a correctable
systematic bias error.

The projects primary accomplishment was the successful application of model-based processing
to signals from micro-cantilever arrays: 40-60 dB improvement vs. the smoother algorithm was
demonstrated. This result was achieved through the development of appropriate mathematical descriptions
for the chemical and mechanical phenomena, and incorporation of these descriptions directly into the
model-based signal processor. A significant challenge was the development of the framework which would
maximize the usefulness of the signal processing algorithms while ensuring the accuracy of the
mathematical description of the chemical-mechanical signal. Experimentally, the difficulty was to identify
and characterize the non-target signals present in the measurement system. In the future, these signals will
limit the ability of the sensor to detect very small quantities of chemicals generated by nuclear processing.
In this project, it became necessary to make use of a model system, mercaptoethanol, which created a large,
reproducible signal that could be readily analyzed with the model-based processor. Further, redundant
cantilevers were examined exclusively: all levers were nominally identically functionalized, and no
control levers were used that did not react to the mercaptoethanol signal. To demonstrate the full utility
of the MBP for chemical sensing, the logical and necessary next steps are (1) verify the physical models
used in this study for a variety of solvents and target molecules (this data has already been obtained as part
of this study) (2) make use of control levers, and (3) extend the experimental library to include low
concentrations of chemical targets of practical interest for sensing nuclear processes.

Outline:

1. Introduction
2. Experimental description
3. Physical chemistry
4. Mechanical modeling
5. Signal processing
6. Summary
7. Appendix A: Derivation of adsorption and desorption parameters

1
1. Introduction: Cantilever Sensor Arrays

Micro-machined cantilevers can function as detection devices when one side is

fabricated to be chemically distinct from the other, as shown in Fig. 1.

functionalized side 0.9 m thick 5 m thick

Lever chip
500 m

Fig. 1: Side view of typical micro-machined cantilever showing cantilever dimensions:

0.9 microns thick, 500 microns long (and 100 microns wide not shown). To function as
a sensor, the top of the cantilever (functionalized side) must be made to be chemically
distinct from the bottom.

Functionalization can be accomplished, for example, by evaporating a thin (~10s

of nm) film of metal such as Au on the top of the chip, then immersing the cantilever chip
in a probe chemical that will bind preferentially to the Au thin film. The lever acts as a
sensor when it is exposed to a second, target chemical that reacts with the probe, since
the reaction causes a free energy change that induces a surface stress at the cantilever
surface. Differential surface stress, , in turn, induces a deflection of the cantilever that
can be measured optically or electronically

The total free energy change at the surface, G, can described as the free energy
due to four contributions: GCANT, the mechanical energy associated with bending the
cantilever, GPOLY, free energy resulting from macromolecular conformational entropy
and nonelectrostatic interactions, GOSM, free energy from osmotic pressure of
counterions near the surface of the cantilever, GELEC, the electrostatic free energy [1].

Eq. 1 G = GCANT + GPOLY + GOSM + GELEC

The free energy change is related to , the stress difference between top and
bottom side of the ith cantilever by:

Eq. 2 i (t ) = Gi (t ) i (t )

where G has units of J/mole and is the sum of all of the contributions to free energy
changes at the surface of the ith cantilever, i(t) is the surface concentration of the
species of interest (in numbers of species per cm2) on the surface of the ith cantilever
[2,3].

2
 The surface concentration of the interacting molecules, i(t), is estimated using
Langmuir kinetics, e.g., the first-order Langmuir kinetics has the following form:

d ()
Eq. 3 = ka c(t ) (1 ) kd max
dt

(t )
Here is , c(t) is the bulk concentration of the target molecule in solution in moles
max
per liter, or [M], ka is the adsorption rate constant in [M]-1 s-1, kd is the desorption rate
constant in cm2 #species-1 s-1, and max is the maximum possible surface concentration of
the species of interest in actual numbers of species per cm2. Please see Appendix A for a
more detailed explanation of Eq. 3.

Finally, differential surface stress in the cantilever induces a deflection on the ith
cantilever, zi(t) , according to Stoneys equation:

4 l
Eq. 4 i (t ) = krect zi (t )
3(1 ) W t

where is Poissons ratio, l, W and t are the rectangular cantilever length, width and
thickness, respectively, krect is the cantilever spring constant.

Cantilevers can function as effective sensors in liquid or gas. In a gas

environment, target molecules bind to an oscillating cantilever and produce resonant
frequency shifts that can be directly linked to the mass of the target. Resonant frequency
shift measurements are difficult to measure in liquid due to damping effects, however.
Cantilever sensors in liquid therefore generally rely on measured DC deflection signals of
non-oscillating levers. This report pertains primarily to cantilever sensors operating in
liquids, although some aspects of the signal processing algorithm may be applicable to
cantilever arrays in gas environments as well.

A limitation on cantilever sensors in liquid is that their signal-to-noise ratio (SNR)

is rather low, often 5:1 or smaller. SNR is expected to be significantly lower in many
operational environments of interest. Further, the reliability of fielded cantilever sensors
ultimately depends on their being incorporated into arrays, which increase system
complexity and can make response interpretation difficult. However, these difficulties
are exactly the type that can be overcome with effective signal extraction techniques such
as the model-based approach, the subject of this invention. In this report, the model-
based signal extraction system for cantilever arrays is described that takes advantage of
the fact that multiple redundant signals are available to the signal processing algorithm.

3
2. Experimental description

In this section, the cantilever detection experimental apparatus will be described,

together with results that demonstrate detection of sulfur-terminated molecules with Au-
coated levers. Non-chemical signals are measured and characterized so that these can be
separated from the chemical signals of interest.

2.1 Experimental description: effects of refractive index differences

An eight-cantilever array chip, obtained from Veeco Instruments, was used for
this experiment. Actual lever deflections were not measured; the goal was to determine
what variability exists in the Scentris optical detection system when fluids of differing
indices of refraction are used in the flow cell. To accomplish this, apparent deflections
were observed as light beams were focused on the rigid base of the chip. The Scentris
system uses super-luminescent diodes (SLD) to generate non-monochromatic light,
which is reflected by the cantilevers into a position-sensitive detector optimized for light
with an average wavelength of about 940 nm (Fig. 2).

Fig. 2: Scentris optical detection system. (Diagram from Scentris Application Notes,
Veeco Instruments.)

The lever chip had 1.5 nm Ti and 20 nm Au evaporated onto the surface.
Immediately prior to being loaded into the flow cell, the chip was cleaned in piranha (1:3,
hydrogen peroxide: sulfuric acid) for 20 minutes, then rinsed in 10 ml water for ~ 10 min.
The chip was then placed directly into the freshly-cleaned flow cell, which was flooded
with water and sealed.

Light was focused on the chip at the base of the cantilever, as shown in Figs. 3
and 4. It appears that the Au surface may have been damaged somewhat by the extended
piranha cleaning (some gold delamination appears to have occurred), but the areas used
for this experiment were smooth and reflective. The actual position of the light spot on
the cantilever is also visible as a white spot in Fig. 3. The positions all of the light spots
used during the experiments are shown by green circles in Fig. 3, and by a green line in

4
the side view of the chip shown in Fig. 4. Note from Fig. 4 that the levers are actually
525 microns thick under the light spots as they were positioned for these measurements.
Due to the fact that there is a damaged area on the surface of the third lever from the top
in Fig. 3, this lever was not used for these experiments.

Lever 1 (SW #)

Lever 2 (SW #)

Lever 1 (SW #)
Lever 3 (SW #)

LE
Lever 4 (SW #)

Lever 5 (SW #)

Lever 6 (SW #)

Lever 7 (SW #)

Lever 8 (SW #)

Fig. 3: Top view of light spot position and cantilevers; levers are numbered according to
software (SW) scheme, reverse of etched numbers on chip.

Fig. 4: Side view of lever chip showing SLD spot position in green.

5
 The Scentris software requires that an effective length for the cantilever be
entered so that deflection signals may be appropriately calibrated. For typical Scentris
experiments, this effective length is measured from an etched mark on the top of the chip
to the center of the light spot, shown schematically in Fig. 3 with the light spot
represented by a yellow dot. (The etched mark is under the red line on the right.) For this
experiment, however, the effective length is not physically meaningful since the light
spot is not being positioned on the flexible part of the cantilever. However, an effective
length of 460 nm was entered into the software so that the magnitude of the calculated
deflection signals might be appropriately interpreted when comparing with actual lever
deflection experiments, where effective lever lengths of about 460 nm are common.

In a first test, ~ 1 ml 5% mercaptoethanol (refractive index: 1.50) in water

(refractive index: 1.33) was manually injected into through the 50 microliter flow cell in
about 10 sec. An apparent deflection occurred, as shown in Fig. 5. The magnitude of the
apparent deflection varied from less than 10 nm for most levers to 20 nm for one lever.

Fig. 5: Effect of injecting 5% mercaptoethanol in water into the flow cell which
previously contained pure water. Nanometer values in the legend are the initial vertical
position of the levers as measured by the Scentris; all subsequent deflections are
relative to this value.

6
In Fig. 6, this experiment is repeated several times; the signals are repeatable and
reversible.

Fig. 6: 5% mercaptoethanol in water causes apparent deflections in optical measurement

system that are repeatable.

7
 In a second test, pure ethanol (refractive index: 1.36) was flowed into the cell.
This caused causes longer-term transients (possibly due to incomplete mixing/exchange
in the flow cell) and somewhat larger apparent deflections for some levers, as shown in
Fig. 7.

Fig. 7: Longer transients and larger apparent deflections (up to ~ 50 nm, downward) for
some levers occur when pure ethanol replaces water in the flow cell.

8
 Finally, it was observed that the optical deflection system is sensitive to ambient
light, even with the opaque hood in place over the test volume. Fig. 8 shows the effects
of turning on and off the fluorescent room lights.

Fig. 8: Small apparent deflections occur when room lights are turned on and off, even
with opaque hood in place over test volume.

9
2.2 Experimental description: deflection due to mercaptoethanol

Experiments were performed with 750 m levers, which were coated by

evaporation with 1.5 nm Ti and 20 nm Au. Levers were cleaned in piranha (3:1 sulfuric
acid: hydrogen peroxide) for ~5 minutes, with gentle agitation, then rinsed in 10 ml water
for ~40 minutes before being loaded directly into the freshly-cleaned flow cell, which
was flooded immediately with purified water and sealed.

The effective length of the levers was 520 microns, as shown in Fig. 9. Lever 6
was broken during chip handling.

Lever 1 (SW #)

Lever 2 (SW #)

Lever 3 (SW #)

Lever 4 (SW #)

Lever 5 (SW #)

Lever 6 (SW #)

Lever 7 (SW #)

Lever 8 (SW #)

Fig. 9: Scentris chip used for mercaptoethanol desorption experiments.

Constant flow experiments in the 50 microliter cell were performed at 50

microliters/min, alternating purified water and 5% mercaptoethanol in water. Results of
the initial steps in the experiment are shown in Fig. 10. Nanometer values in the legend
are the initial vertical position of the levers as measured by the Scentris; all subsequent
deflections are relative to this value for each lever. Note that levers with similar initial
deflection values tend to drift at similar rates throughout the experiment. The reasons for
this behavior are not yet clear; it could be related to the mechanical properties of the
levers, or it could be an artifact of the optical detection system.

10
Fig. 10: Initial stages of mercaptoethanol desorption experiment.

An initial positive (downward, away from the Au side) deflection occurs

following the first injection of 5% 2-mercaptoethanol in water. This behavior is
consistent with an increase in compressive stress in the cantilevers. This initial positive
deflection has been observed repeatedly in previous experiments with mercaptoethanol
and Au-coated levers. Subsequent binding of mercaptoethanol to the Au surface induces
tensile stresses in the lever, manifested by a negative (upward, toward the Au) deflection
of the levers.

When water is flowed through the cell, additional compressive stress is developed
in the cantilevers; this is shown in two events at ~8200 sec and 12000 sec in Fig. 10. In
Fig. 11, the water/5% mercaptoethanol exchange is repeated twice more, and behavior is
consistent: levers bend downward (decrease in compressive stress) when

11
mercaptoethanol is introduced, then relax upward (increase in compressive stress) when
pure water flows into the cell. Note that the levers never fully relax upward to their
initial positions, indicating that significant numbers of mercaptoethanol molecules remain
bound to the Au lever surfaces even after rinsing for long periods.

Fig. 11: Final stages of mercaptoethanol desorption experiment.

12
2.3 Experimental description: effects of fluid flow

750 micron long Au levers (1.5 nm Ti, 20 nm Au, by evaporation), on the same
chip used for the mercaptoethanol experiments previously described, were cleaned in
piranha (3:1 sulfuric acid: hydrogen peroxide) for 5 minutes, then rinsed in 10 ml water
for 55 minutes. They were loaded directly into freshly-cleaned flow cell, which was
flooded with water. The SLD spot positions were set to maintain the levers effective
length near 520 m. As can be seen in Fig. 12, some pitting of the Au surface may have
occurred since the previous experiment, possibly due to the piranha cleaning process.

Fig. 12: Levers used for ethanol desorption and flow rate experiments

Constant flow experiments were initially performed in the 50 microliter flow cell.
Water and 5% ethanol in water were exchanged several times. Relative to previously-
observed deflections induced by 5% mercaptoethanol in water, the magnitude of these
deflections were small, as shown in Fig. 13. The behavior of the levers following each
5% ethanol injection was consistent and reproducible. These experiments demonstrate
that the large deflections observed with the mercaptoethanol in previous experiments
were likely caused mainly by molecules adsorbing or desorbing from the surface, as
opposed to being caused by changes in water-Au surface interactions mediated by
ethanol.

13
 1000
lever 3 (SW #): 498 nm
Start 5% ethanol in water, 50 ul/min lever 4 (SW #): -1335 nm
lever 5 (SW #): -322 nm
lever 7 (SW #): -787 nm
800 Water at 50 microliters/min lever 8 (SW #): 3524 nm
Avg

5% ethanol in water at 50 ul/min

600

400
Deflection (nm)

200

-200
Water at 50 ul/min

5% ethanol in water at 50 ul/min

-400
Start water flow at 50 microliters/min Water at 50 ul/min

3
0 5 10 15 20x10
Time (s)

Fig. 13: Constant flow experiments with water and 5% ethanol in water.

It is interesting to note in Fig. 13 that there is an initial positive deflection of the

levers by ~175 nm when flow is initiated at 50 microliters/min. This behavior is
important to understand so that flow-induced deflection is not mistaken for deflection
induced by chemical changes at the lever surface. Subsequent experiments were
therefore performed to learn more about these flow-induced deflections. In Figs. 14 and
15, deflections are shown as the flow rate is slowed from 50 to 2 microliters/min.

14
 40 0

300
t=21300sec : Change flow rate to 20 microliters/min

t=21620sec : Change flow rate to 10 microliters/min

t=21936sec : Change flow rate to 5 microliters/min

200
t=22332sec : Change flow rate to 2 microliters/min
Deflection (nm)

10 0

lever 3 (SW #): 498 nm

0 lever 4 (SW #): -1335 nm
lever 5 (SW #): -322 nm
lever 7 (SW #): -787 nm
lever 8 (SW #): 3524 nm
Avg
Water at 50 ul/min

-100
3
20 25 30 35 40x10
Time (s)

Fig. 14: Effect of reducing water flow rate from 50 to 2 microliters/min.

The levers drift downward, then back up after the flow rate is stabilized at 2
microliters/min for an extended period. Fig. 15 shows an expanded view of the average
lever deflection as the flow rate is changed from 20 to 2 microliters/min.

15
 t=21300sec : Change flow rate to 20 ul/min from 50 ul/min

38
t=21620sec : Change flow rate to 10 microliters/min

36 t=21936sec : Change flow rate to 5 microliters/min

34 t=22332sec : Change flow rate to 2 microliters/min

Deflection (nm)

32

30

28

lever 3 (SW #): 498 nm

26 lever 4 (SW #): -1335 nm
lever 5 (SW #): -322 nm
lever 7 (SW #): -787 nm
lever 8 (SW #): 3524 nm
24 Avg

3
21.5 22.0 22.5x10
Time (s)

Fig. 15: Average lever deflection as water flow rate is decreased to 2 microliters/min
from 20 microliters/min.

In Fig. 16, some deviations occur in the downward drift trend as the flow is
reduced, but these are relatively small. Fig. 16 shows the effect of increasing flow rate
from 2 to 200 microliters/min, then back to 2 microliters/min. It is difficult to generalize
about the direction of deflection induced by increased or decreased flow. Figs 14 and 15
show downward deflection as flow rate is decreased, while Fig 16 shows downward
deflection as flow rate is increased. In general, however, larger flow rates do seem to
cause larger deflections, and when the flow rate is stabilized for extended periods at low
values (<5 microliters/min), the levers tend to drift upward for this experimental set.

16
Fig. 16: Lever deflections as water flow rates are changed.

Fast injection conditions are of particular interest for detection experiments.

Ideally, chemicals can be introduced rapidly into the fluid cell so that the resulting
behavior of the levers can be characterized starting at a single time (as opposed to the
case when chemicals are slowly flowed into the cell over minutes or hours). Also, since
high flow rates themselves can cause lever deflections, the ideal chemical detection
experiment would occur under low- or no-flow conditions. Data in Fig. 17 show the
results of changing flow rates briefly from very low levels (2 microliters/min), briefly to
high levels for injection (100 or 200 microliters/min), then back to low levels for
detection.

17
Fig. 17: Effect of increasing water flow rates briefly for simulated injection of a
chemical of interest into the flow cell; the initial flow rate is 2 microliters/min.

As can be seen in Fig. 17, increasing flow rates to 100 or 200 microliters/min
from 2 microliters/min causes a large, immediate upward deflection. More problematic
for detection experiments, is that the fast changes can cause extended relaxation periods
(1000s of seconds) during which the levers drift significantly.

18
2.4 Experimental description: temperature and chemical effects

To calibrate the deflection of Scentris cantilevers, it is helpful to use a

temperature stimulus which can be measured independently. As will be described in
Section 4, the spring constant of the cantilever can be calculated based upon the known
dimensions of the cantilever and the observed deflection response of the cantilever to an
independently measured temperature increase. As can be seen in Fig. 18, critical
dimensions of the cantilever can be determined by scanning electron microscopy (SEM).

Fig. 18: SEM image of cantilevers showing sensitive measurement of the cantilever
thickness.

Here, we describe a calibration experiment in which 500-micron long, Au-coated

levers were exposed to a temperature pulse, then a chemical stimulus in the form of a 1%
(by volume) concentration of mercaptoethanol in pure water, and finally a second
temperature pulse.

Levers were coated with 20 nm of evaporated Au, and cleaned as described

previously with piranha immediately before being loaded into the water-filled 50
microliter flow cell. Except for brief pauses to exchange syringes in the syringe pump,
the entire experiment was performed at a constant flow rate of 20 microliter/min. The
results of the experiment are shown in Fig. 19.

19
 1200
Deflection (nm)
Temperature Change (K) 10
1000
mercaptoethanol

Temperature Change (K)

800 injection at 20 ul/min 8
Deflection (nm)

600
6
water injection at 20 ul/min
400
4
200

2
0

-200 0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Time(sec)

Fig. 19: Deflection of cantilever due to temperature and chemical stimuli.

The results of this experiment will be discussed in detail in Sections 3-5, where
the deflections are modeled and quantitatively analyzed to extract chemical parameters,
such as the absorption constant, that are necessary for detection applications.

In order to ensure that experimental results are not being affected by uncontrolled
variables, and that modeling assumptions are realistic, a large data set was obtained with
various concentrations, target molecules and solvents. The experimental parameters of
these experiments are summarized in Table 1, and representative results of the
experimental data set are shown in Figs. 20 and 21.

Table 1: Number of experiments for several experimental conditions:

Solvent: Water Ethanol Ethanol
Solute: Mercaptoethanol Mercaptoethanol Hexadecane
Concentration (M): 1.E-04 0 0 2
1.E-03 0 6 3
1.E-02 1 1 0
1.E-01 0 3 0
1.E+00 3 2 0

20
 1000
040926_Au_EthanolHex_001: 1 mM hexadecane thiol in ethanol
8

500
6

Delta temp (C)

Deflection (nm)

4
0 Flow pure ethanol

2
1 mM hexadecane thiol in ethanol
-500

-1000
3
60 70 80x10
Time (s)

Fig. 20: Experimental results for 1 mM hexadecane in ethanol.

1500
10

040916_Au_Ethanol_000: 10 mM mercaptoethanol in ethanol

1000 8

6
500
Delta temp (C)
Deflection (nm)

0
2
flow pure ethanol
Flow 10 mM mercaptoethanol
0
-500

-2

-1000
3
50 55 60 65 70 75 80x10
Time (s)

Fig. 21: Experimental results for 10 mM mercaptoethanol in ethanol.

21
3. Physical chemistry

In this section, we describe how we predict changes in surface stress as a function

of surface loading. The overall analysis in the modeling begins with Stoneys equation
(Eq. 4), which states that the deflection of the cantilever is directly proportional to the
difference in surface stress on the cantilever surface. The physical system is shown
below in Fig. 22.

Fig. 22: Flow over a cantilever with polymer, green species, depositing on the surface of
a cantilever. The deposited polymer self-assembles and the assembled species induces a
lateral stress that causes a stress imbalance between the top and bottom surfaces, i.e., +
and -.

This stress differential constitutes the signal. We follow the approach as

described by Lavrik et al. [3] to relate the surface stress difference to the surface coverage
and the free energy of adsorption (Eq. 2). What remains is to develop physical models
for (t) and G .

To estimate the rate of change of surface concentration, (t), we employ a

modified form of Eq. 3 by re-defining a rate of desorption, kd' , where kd' = kd max . With
this new definition, Eq. 3 becomes:

22
Eq. 5
d ()
= ka c(t ) (1 ) kd'
dt

All terms are as described above in Eq. 3 with the exception of kd' . Note that Langmuir
kinetic equations have been presented in various forms in the open literature. The above
form, Eq. 5, was chosen to enable a one to one comparison between rate constants
predicted by the MBP presented here and those reported in the literature [4,5,6]. Eq. 3
above, however, is based on the rigorous analysis presented in Appendix A.

In solving the Langmuir kinetics, we hold the concentration as a constant in Eq. 3

and then build in the dynamic nature of the signal through the solution to a decoupled
differential equation for the dynamics of the bulk phase signal concentration (see Eq. 7).
This is justified since the time scales between surface adsorption/desorption and bulk-
signal concentration are sufficiently disparate.

Once the input signal, c0, is turned off we again fit the desorption process with an
nth order Langmuir desorption mode:

Eq. 6
d ()
= kd ( )
n

dt

The fit to the actual desorption data revealed a third-order dependence on the
dimensionless surface concentration, or n = 3. Desorption processes are typically
modeled with a first or second order model to extract kd; however, in our case, the best fit
was with a third order model; consequently, the third order dependence on dimensionless
surface concentration was not expected. This fact that the best fit for the desorption
behavior was with n>1 suggests that there may be a distribution of activation energies
associated with the desorption process (see Appendix A).

To account for the dynamics of the system, more specifically the transient nature
of the system owing to fluid flow rates, we developed a dynamic estimate of the bulk
species concentration based on a modified stirred tank reactor system, see Fig. 23 below.
This mathematical description of the flow cell is a useful approximation for observing
time-dependent fluidic behavior, but experimental evidence suggests that at the low flow
rates applied in the 50 ul flow cell experiments described here (50 microliters/min and
smaller), flow is in fact relatively laminar. Levers on the outside of the array, for
example, were consistently observed to react to chemical signals a few seconds before the
inner levers. If the system acted as a true stirred tank, all levers would react to chemical
stimuli simultaneously.

23
Fig. 23: Model flow cell for cantilever array. Qin and Qout are the volumetric flow rates
in and out of the cantilever array chamber. ci and c are the input concentration and the
concentration within the volume for perfect mixing.

The resulting equation we used to predict the bulk, target species concentration
was the following form:

Eq. 7

{( ) ( ) }
c(t ) = c0 1 e (t tON ) / H ( t tON ) 1 e ( t tON ) / H ( t tOFF )

Here c(t ) is the time dependent concentration of our target species, ci, is the input signal
or concentration, tON is the time at which the signal is turned on and tOFF is the time when
V
the signal is turned off and , = , is the residence time of the fluid in the cantilever
Q
array volume, V. H(t) here is the Heaviside step function; i.e. the function that is zero
when the argument is negative and one when it is positive. For the given system and
using half of the chamber volume, 75s , and the resulting response or concentration in
the chamber, for the gross approximation of perfect mixing, is shown below in Fig. 24.

24
Fig. 24: Bulk concentration in cantilever array flow chamber. Given that the cantilever
array is at the center of the chamber, predictions were made for the entire chamber volume,
blue curve, and for half of the chamber volume, brown curve. The broad black line
represents the time on and the time off of the signal.

All of these features were built into the model based signal processor (Section 5),
and the signal based processor/estimator was used to fit these coupled equations predict
the optimal values for ka, kd and max. With these parameters specified, we calculate the
k
free energy for adsorption, G = RT ln a [4,6]. For any target species, we can then
kd
fit to get the appropriate adsorption parameters. Furthermore, we also built in the
cantilever response to thermal effects, or coefficient of thermal expansion effects, as will
be described in Section 4.

4. Mechanical modeling

In this section we describe the modeling of micromachined cantilevers using finite

element methods and closed form analytical expressions. The results are applied to
calibrate the cantilevers based on their response to known temperature stimuli.

4.1 Finite element model of a complete array

A finite element model of an eight-cantilever array was created to conduct

simulations. The wafer, cantilevers, and anchors are represented by eight-node solid
elements. The gold coating on the cantilever surfaces are represented by four-node
membrane (shell with one integration point) elements. The boundary conditions are fixed
(no displacement, no rotation) on the bottom surface of the chip. This will account for
compliance at the root of the cantilevers.

25
Using this model, various loading conditions can be applied, e.g., fluid interaction, heat
transfer, applied surface traction. Dynamic behavior can then be simulated using the 3D
nonlinear, implicit finite element code, NIKE3D. The model is sufficiently generalized to
account for variations in geometry, material properties, and initial stress.

Fig. 25: Finite element mesh of an IBM Zurich Research A1: Water, cantilevers, and
anchor.

4.2 Finite element model of a single cantilever

A finite element model was created to focus on the behavior of a single cantilever.
As in the full array model, the cantilever, anchor, and a portion of the wafer are
represented by eight-node solid elements. The gold coating on the cantilever is
represented by four-node plate elements. In this case, only membrane behavior is
captured by using one integration point with the shell elements.

Two mechanisms which cause cantilever deflection were studied with the model.
The first is temperature-dependent thermal expansion. The cantilever deflects in the same
fashion as bi-metallic beams. As the temperature changes, the silicon and gold expand at

26
different rates (as described by the thermal expansion coefficient). This causes a stress in
the beam resulting in a deflected equilibrium shape.

The second mechanism is a chemically induced stress which also causes the
cantilever to deflect to maintain equilibrium. Chemical adsorption/desorption in the gold
layer causes a change in the free energy of the gold layer. This results in an applied
traction to the top surface of the silicon cantilever. To maintain static equilibrium the
cantilever must change its shape. Due to the beam-like geometry of the cantilever, the
deflection is predominately a displacement in the direction normal to the cantilevers
surface (z direction in Fig. 26).

This effect cannot be captured directly in the current version of NIKE3D,

However, it can be included by an equivalent thermal stress in fictitious shell elements.
Those shell elements are given properties so that they do not contribute to the mass and
stiffness of the cantilever and that plane stress arises from contributions from thermal
expansion but not bending, i.e., the thermal expansion coefficient, , is very large relative
to the Youngs modulus. From the plane stress equation,

2 (xx + yy )
E E
Eq. 8 xx = T
1 1+

if is very large, then xx is dominated by the second term on the right hand side. In
these equations, E and are the Youngs modulus and Poissons ratio for Si,
respectively. T is the change in temperature. xx and yy are the in-plane strains.

Thermo-elastic material properties (Material Type 4 in NIKE3D) are taken to be

for silicon: E = 150e9 kg/m3, = 0.2, and = 4.2e-6 / K and for gold: E = 80e9 kg/m3,
= 0.42, and = 14.2e-6 / K.

Fig. 26: Finite element mesh of a single cantilever.

27
Fig. 27: Example of the deflection shape of a single cantilever under loading. Contours
are of displacement in the normal (z) direction.

4.3 Analytical relationships

4.3.1 Deflection of a bi-metallic cantilever due to a change in temperature. From Roark

and Young, Formulas for Stress and Strain, Fifth Ed., Sect. 7.2[7],

Eq. 9
ztip 2 +
= 3L Au 2 Si
T Si k1

where ztip is the normal (z-direction) deflection of the free end (tip) of the cantilever, T
is the temperature difference, is the difference in thermal expansion coefficients
(Si Au), L is the cantilever length, is the thickness of the Si or Au, and k1 is

Au 2 E (1 Si ) Au 3 E Si (1 Au) Si
k1 = 4 + 6 + 4( Au ) + Au ( ) + ( )
Si Si E Si (1 Au ) Si E Au (1 Si ) Au

where E is Youngs modulus for Si or Au and is Poissons ratio for Si or Au. Note the
Poisson effect has been included in k1.

4.3.2 Deflection of a cantilever beam from a resultant surface stress. A modified version
of Stoneys equation (applies free-free beam), which accounts for the fixed end boundary
condition, comes from, Sader [8].

28
Eq. 10
1 1 2 E
= z
K 3 L 1

where
1 1 1 1 1 1 2
K = X + 2X + B + + 2 2 + 2 + B
2
1 2 12 1 2 1 2

and

[
1, 2 = 2 3 5(1 ) m 10(1 )(2 3 ) ]
In these equations, is the surface stress (J/m^2 or N/m), is the thickness of
the Si cantilever, L is the length of the cantilever, E is Youngs modulus of Si, is
Poissons ratio for Si, z is the normal (z-direction) deflection of the free end (tip) of the
cantilever, X = x/L = 1 is the normalized length, and B = w/L is the normalized width of
the cantilever. Note that , surface stress, is related to the surface stress energy per unit
area for this cantilever geometry. (For details, see W. Haiss, [9])

K is the correction factor for Stoneys equation. For the case, L = 500 microns, w = 100
microns, and t = 1 micron, K = 1.04

4.3.3 Spring stiffness. A cantilever beam can be characterized by a spring stiffness

analogous to a one-degree-of-freedom mass/spring oscillator. The deflected shape
corresponds to the first mode of vibration,

Ew
3

Eq. 11 k=
4 L

where E is Youngs modulus, w is the width, is the thickness, and L is the length. The
moment of inertia of the cantilever with respect to the normal direction (often seen in
literature) is

1
Eq. 12 I= w 3
12

4.3.4 Deflected shape of the cantilever. The normal displacement at any point along the
cantilever is

Eq. 13 z ( x) =
f
(
x 2 x 2 + 6 L2 4 Lx )
2 Ew 3

29
where x is the position along the cantilever (x=0 is the fixed end. x=L is the free end). It
is assumed that the cantilever is subject to a uniform loading, f (force/length).

The ratio of displacement at any point on the cantilever with respect to tip
displacement is then

z(x) x 2 (x 2 + 6L2 4Lx)

Eq. 14 =
ztip (L) 3L4

This relationship can be used when the measurement point is not located at the tip of the
cantilever. For example, given a 500 micron cantilever, the displacement at a point 355
microns from the free end is 61.57% of the tip displacement.

4.4 Finite element calculations vs. analytical relationships

The following section uses the material properties values: Silicon: E = 150e9 Pa,
= 0.2, = 1.0e-6 m, L = 500.e-6 m, = 4.2e-6 / C. Gold: E = 80e9 Pa, = 0.42, =
20e-9 m, L = 500.e-6 m, = 14.2e-6 / C.

Consider tip deflection due to temperature difference: the analytical expression in

Section 4.3.1 yields z / T = 103 nm/K. From the finite element model, z / T = 111
nm/K. The analytical relationship is one dimensional. It does not take into account
shortening in the axial direction. The finite element model does capture this effect. Figure
28 shows the corresponding amount of shortening in the axial direction for a 10 degree
temperature increase.

30
 -7
x 10
1.4

1.2

Vertical Displacement (m)

0.8

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Temperature Increase (K)

Fig. 28: Vertical displacement vs. temperature increase. Red: finite element model. Blue:
analytical relation.

-8
x 10
2.5

2
Axial Displacement (m)

1.5

0.5

0
0 1 2 3 4 5 6 7 8 9 10
Temperature Increase (K)

Fig. 29: Axial displacement vs. temperature increase from finite element model. The
axial displacement is two orders of magnitude smaller than the vertical displacement.

31
4.5 Calibration of cantilevers using temperature-deflection-determined spring constant

Material, geometric, and manufacturing (i.e., residual stress) variations in the

cantilevers cause variations in mechanical response. Experimental measurements show
that individual cantilevers respond to a given temperature change with varying
displacement. We have attempted to characterize these variations through the spring
constant (implemented for the purposes of our calculations in the Youngs modulus) of
each cantilever.

Experimentally, a known temperature increase in measured along with the

corresponding cantilever displacement. Using the relationship in Section 4.3.1, a
corresponding spring constant, k, and a Youngs modulus, E, can be calculated. The
ideal value of dz / dT is 104 nm/K which corresponds k = 30 nN/m.

3000

E = 50 GPa: dz/dT = 261 nm/K

2500

2000
Vertical Deflection (nm)

1500 E = 100 GPa: dz/dT = 140 nm/K

1000
E =150 GPa: dz/dT = 96 nm/K

500 E = 200 GPa: dz/dT = 73 nm/K

0
0 1 2 3 4 5 6 7 8 9 10
Temperature Difference (K)

Fig. 30: Vertical displacement vs. temperature increase from four experimental
measurements. For each dz / dT, an experimentally determined E is calculated.

An example of temperature vs. deflection data for eight cantilevers is shown in

Fig. 31. A least-squares-fit line was used to estimate dz/dT for each lever. Then using the
relationships in Sections 4.3.1 and 4.3.3, the corresponding spring constants were
determined for each lever. The values ranged from 12.2 mN/m to 13.4 mN/m.

32
 8 lever example: hexadecane and ethanol
100

50

0
Vertical Displacement (nm)

-50

-100

-150

-200

-250
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Temperature Increase (C)

Fig. 31: Vertical displacement vs. temperature increase from a single temperature pulse.
Displacement response of eight cantilevers. Spring stiffnesses varied from 12.2 mN/m to
13.4 mN/m.

As an example, consider the data shown in Fig. 19. Here, we can fit to get the
appropriate adsorption parameters. Furthermore, we also built in the cantilever response
to thermal effects, or coefficient of thermal expansion effects. To estimate the
contribution of CTE, we use Eq. 9 and find that z/T approximately equals 98 nm/K.

The deflection, z, is a function of material properties and temperature of the

system. This estimate was based on literature values for material properties, and later in
the project, the estimate was refined through temperature calibration experiments for each
cantilever in the array. With all of these features built into the model based signal
processing software, with fit parameters, we now have a capability to find a signal, with
the appropriate features that our target species would produce, that is buried in non-signal
or noise.

5. Model-based signal processing of cantilever signals

This section describes model-based processing, and how it is applied to achieve an

enhanced signal for the cantilever array in a chemical detection problem.

5.1 Model-based signal processing approach

The model-based approach is incorporating mathematical models of both

physical phenomenology and the measurement process (including noise) into the

33
processor to extract the desired information. This approach provides a mechanism to
incorporate knowledge of the underlying physics or dynamics in the form of
mathematical process models (along with measurement system models and
accompanying noise as well, as model uncertainties) directly into the resulting processor.
In this way, the model-based processor enables the interpretation of results directly in
terms of the problem physics. The model-based processor is actually a modeler's tool
enabling the incorporation of any a priori information about the problem to extract the
desired information.

The fidelity of the model incorporated into the processor determines the
complexity of the model-based processor with the ultimate goal of increasing the inherent
signal-to-noise ratio (SNR). These models can range from simple, implicit, non-physical
representations of the measurement data such as the Fourier or wavelet transforms to
parametric black-box models used for data prediction, to lumped mathematical
representations characterized by ordinary differential equations, to distributed
representations characterized by partial differential equation models to capture the
underlying physics of the process under investigation. The dominating factor of which
model is the most appropriate is usually determined by how severe the measurements are
contaminated with noise and the underlying uncertainties. If the SNR of the
measurements is high, then simple non-physical techniques can be used to extract the
desired information; however, for low SNR measurements more and more of the physics
and instrumentation must be incorporated for the extraction.

This paper is concerned with the novel development and application of model-
based signal processing to micro-machined cantilever sensor arrays. We will first
functionally describe the generic approach in terms of the cantilever system and then we
will develop the cantilever and associated mathematics more specifically eventually
developing a first-cut design based on these models. The overall approach captures all
possible models in terms of nonlinear differential equation models with the one described
specifically an embellishment of the model-based approach.

We briefly describe the generic model-based approach to the cantilever array problem as
depicted in Fig. 32. The movement of cantilevers is influenced by three major factors:
chemical binding, temperature effects and fluidic forces. We will assume in this
formulation that the temperature dynamics are relatively slow. Therefore, the changes in
temperature can simply be modeled by a time-varying polynomial. Fluidic disturbances
act more quickly, and can be described by the linear Stokes equations of motion. As a
result, one can also solve the fluids problem separately. The resulting predicted stress
can be superimposed over stresses due to temperature and chemical effects. Chemical
reactions at the cantilever surface act over tens of seconds to many minutes, and represent
the signal of interest. The entire process is contaminated with both measurement and
process noise, which decreases the overall signal-to-noise ratio (SNR) and makes this
invention so important. The poor signal levels can be enhanced by incorporating
mathematical models of the underlying physics and measurements including the noise
statistics into the processor to produce an enhanced signal (deflection) with a much
higher SNR.

34
 PHYSICS
Thermodynamic,
fluidic equations

NOISE Model-based Enhanced signal

Random and processor
deterministic

MEASUREMENT
Deflection

Fig. 32: Generic Model-Based Processing for Cantilever Sensor Arrays.

Mathematically, the time-dependent differential surface stress, (t), on each

cantilever is a function of the surface concentration given by (t ) , which is in turn
parameterized by the maximum possible surface concentration, max , and the initial
target molecule concentration c0 and time-dependent target concentration c(t) in solution.
The time rate of change of the surface concentration, (t ) , can be expressed as

Eq. 15 (t ) = f [(t ),c0 , c(t ),ka ,kd , max ] + w(t ) ,

where w(t) is time-dependent process noise, assumed to be Gaussian, and all other terms
are as defined above.

The change in a cantilevers surface free energy, G(t), is related to the surface
concentration and the target molecule concentration in solution by

Eq. 16 G (t ) = f G [(t ), C (t ), (t ), ka ,kd ]

where (t ) is the time-dependent temperature. The differential surface stress is a function

of the change in free energy and surface concentration by

Eq. 17 (t ) = f [G (t ), (t )]

With this underlying physical phenomenology we can now express the deflection
at the ith cantilever as a function of both surface concentration and induced surface stress
as

Eq. 18 zi (t ) = f z [ i (t ), i (t )]

35
Here we assume that the differential surface stress and as the surface concentration are
variable. The actual measurement is somewhat more complicated due to instrumentation
noise v(t ) as well as disturbances (t ) such as systematic errors, temperature drifts and
just trends. Therefore, we define the cantilever array deflection measurement at the ith
lever as

Eq. 19 di (t ) = f di ( i , zi , i ) + vi (t )

with i a gain or attenuation assumed constant. In this characterization, we also assume

that the additive noise processes are zero-mean and Gaussian. This representation can be
further parameterized in terms of an approximate Gauss-Markov model with nonlinear
dynamics and algebraic equations. We will discuss this subsequently. This completes the
generic description of the system cast into the model-based framework, next we define
some more specific representations based on the underlying phenomenology.

36
5.2 Model-based processor design for cantilever arrays

In this section we develop the model-based approach first for the generic model
sets and then for the specific embellishment discussed previously. We start the section
with the development of an approximate Gauss-Markov model, which can be used to
capture the cantilever signal enhancement problem and then specifically apply it to the
micro-cantilever array sensor system.

If we have nonlinear dynamics (differential equations) describing the system

under investigation, then an approximate representation of the deterministic process and
associated measurement is easily captured in state-space form. State-space is simply
converting an nth-order set of coupled differential equations into an equivalent set of n-
first order differential equations. With this accomplished, we obtain the general nonlinear
vector functional relations defined by the process and measurement models,

x& (t ) = a [ x, u; ] [state]

Eq. 20
y (t ) = c [ x, u; ] [measurement]

where x is the N x -dimensional state vector; y is the N y -dimensional measurement. If

these processes are contaminated by additive zero-mean, Gaussian noise processes then
the approximate Gauss-Markov model evolves as

x& (t ) = a [ x, u; ] +w(t) [process (state)]

Eq. 21
y (t ) = c [ x, u; ] +v(t) [measurement]

Here a[], c[] are the N x -dimensional vector process function and N y -dimensional
measurement functions, respectively, for the process noise given by w ~ N (0, Rww ) and
the corresponding measurement noise as v ~ N (0, Rvv ) . N (, ) is the Gaussian
distribution specified by mean and covariance. With this representation in mind, we can
now define the cantilever signal enhancement problem as:

GIVEN a set of noisy displacement measurements, { y (t )} with known inputs, {u (t )} and

parameters, { } specified by the approximate Gauss-Markov model of Eq. 21, FIND the
best (minimum error variance) estimate of the displacement and surface concentrations,
y (t ), x(t)
, respectively.

37
The solution to this problem can be derived in a wide variety of approaches. Here, we
will use the common Bayesian approach[10]. We summarize the algorithm as:

MODEL-BASED PROCESSOR ALGORITHM

x& (t | t 1) = a [ x , u; ] [State Prediction]

y (t | t 1) = c [ x , u; ] [Measurement Prediction]
Eq. 22

(t ) = y (t ) y (t | t 1) [Innovation or Residual]

x (t | t ) = x (t | t 1) + K (t ) (t ) [Correction]

where K (t ) is the gain of the processor, which must be calculated from the underlying
process statistics (see Candy, Signal Processing: The Model Based Approach
(1986)[10]) for details.

For the micro-cantilever array problem, we must convert the physical relations,
that is, Eqs. 2-4 into the state-space form above. We chose to solve the differential
equation and incorporate the resulting relations into the measurement model. For the state
we modeled the free energy as a piecewise constant function, converted it to discrete-time
using the first difference approximation and excite it with zero-mean, white Gaussian
(process) noise creating a random walk model for this parameter. Therefore, we start with
defining the state vector as, x := G and the measurement, y := d , then we obtain the
following relations. From Eq. 22, the surface concentration relation becomes

Eq. 23 x(t ) = x(t 1) a[ x; ] [Process Model]

Now the measurement model is tricky. We must first solve for the physical variables to
obtain the generic form of Eq. 21; therefore, we have from Eqs. 2-4 that

Eq. 24 yl (t ) = d l (t ) + (t ) c[ x, u; ] [Measurement Model]

The dynamic surface concentration, (t) for adsorption-desorption and pure desorption
processes that were used have the following form:

38
Eq. 25

0
t < tON


(t ) =
(
1 e( kaC (t )+ kd )( t tON )
) for tON t tOFF


1 t > tOFF
2kd ( t tOFF )

c(t ) k
Here = ; K EQ = a , and tON denotes when the signal is introduced and
c(t ) + K EQ kd

tOff denotes when the signal is turned off.

As described above, the differential surface stress on the cantilever is a function of the
surface concentration and free energy

Eq. 26 (t ) = (t )G (t ) M 1 = (t ) x(t )

Therefore the deflection of the lth -lever is weighted by the Stoney equation with a
specific modulus of the individual lever, that is,

3L2
Eq. 27 d l (t ) = l (t ) for l := (1 )
El

where l is the lth Stoney equation with cantilever modulus, El , with lever length and
thickness, L, , respectively and is Poissons ratio for silicon. Thus, we obtain the
measurement equation at the lth -cantilever as

Eq. 28 yl (t ) = l (t ) x(t ) + (t ) for l (t ) := l (t )

and (t ) is the known temperature profile (low-pass filtered). Finally, assuming that both
noise sources are Gaussian random processes (as before), then the result is a time-varying
Gauss-Markov (not approximate due to linearity) multi-channel cantilever model defined
by

39
 x(t ) = x(t 1) + w(t 1)

Eq. 29
yl (t ) = l (t ) x(t ) + (t ) + vl (t )

for w ~ N (0, Rww ) ; v ~ N (0, Rvv ) . Then as before, we can develop the model-based
processor based on this Gauss-Markov model. First, we define the signal enhancement
problem in terms of the cantilever models as:

GIVEN a set of noisy N y -vector displacement measurements, {y (t )} with known

N u -vector inputs, {u(t )} and parameters, { k } specified by the Gauss-Markov
model of Eq. 27, FIND the best (minimum error variance) estimate of the
displacement and N x -vector surface concentrations, y (t | t 1), x (t | t ) ,
respectively.

The model-based algorithm to solve this problem using the specified models is:

Eq. 30

CANTILEVER ARRAY MODEL-BASED PROCESSOR ALGORITHM

x& (t | t 1) = x (t 1| t 1) [Surface Concentration Prediction]

y l (t | t 1) = l (t )x (t | t 1) + (t ) [Displacement Prediction]

l (t ) = y l (t ) y l (t | t 1) [Innovation or Residual]

x (t | t ) = x (t | t 1) + k (t )(t ) [Surface Concentration Correction]

This completes the development of the MBP algorithm for Cantilever Sensor
Arrays. Note that once this framework is developed, it is quite easy to define other
problems of high interest (e.g. detection problems [11]).

5.3 Model-based processor performance evaluation

In this section we discuss the performance of the model-based processor (MBP)

for signal enhancement of an L-element cantilever sensor array. The basic approach we

40
take for MBP performance evaluation is illustrated in Fig. 33. After obtaining the average
parameters by performing the parameter estimation, a Gauss-Markov simulation was
designed to generate synthesized cantilever deflection measurements using the model
discussed in the previous section. Once synthesized at a particular signal-to-noise ratio
(SNR), the processors were applied to the data and their performance analyzed based on
the truth deflections generated by a noise-free simulation. Performance metrics are
applied to evaluate and compare performance. We discuss the various steps in this
procedure.

Parameter Estimation

The basic approach we use is to first parameterize the cantilever array model by
performing parameter estimation on the raw deflection measurements to extract the
critical absorption, desorption and maximum concentration, that is,
l = {ka (l), kd (l), G (l)} ; l = 1,L , L . The parameter estimator we employed was a
nonlinear least-squares criterion using the Nelder-Meade polytope search algorithm [12].
This algorithm is based on minimizing

Nt
Eq. 31 min J ( ) = l2 (t ; ) for (t ; ) := yl (t ) y l (t ; ) ,
l
t =1

where the estimated or filtered cantilever measurement at the lth -lever is given by

Eq. 32 y l (t ; ) = dl (t ; ) + (t ) .

Once these parameters are extracted from the data, they are averaged to give,
{ }
= ka , kd , max . These are the parameters that are used in the Gauss-Markov simulation
model.

We ran this parameter estimator on raw experimental deflection data and

estimated the parameters for each lever. The results are shown in Fig. 34 where we see
the fitted deflection responses compared to the measured. It is clear that the extracted
parameters reasonably predict the filtered cantilever response of Eq. 31. Next we
investigate the development of the simulator.

We used a nonlinear least-squares algorithm in MATLAB to estimate the parameters.

41
 SIMULATION
Cantilever
Measurements

Parameter
Estimation

Cantilever
MODEL

Synthesized
Deflection
Measurements

Smoothing Model-Based
Processor Processor

Error
Statistics

Performance
Analysis
ANALYSIS

Fig. 33: Model-Based Processor Performance Evaluation: Simulation, Processing and

Analysis.

1000 1000
Cantilever No. 1 Cantilever No. 2
800 800
Deflection (nm)

600 600

400
400

200
200

0
0

0 1000 2000 3000 4000 5000

0 1000 2000 3000 4000 5000
Time (sec) Time (sec)
1200 1200
Cantilever No. 4
1000
Cantilever No. 3
1000
800
Deflection (nm)

800
600
600
400

200 400

0 200

0 1000 2000 3000 4000 5000

0 1000 2000 3000 4000 5000
Time (sec)
Time (sec)

1000 1000
Cantilever No. 6
Cantilever No. 5
800 800

600
Deflection (nm)

600
400
400
200

200
0

0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000
Time (sec) Time (sec)

Fig. 34: Parameter Estimator Results: Predicted response (solid line) compared to raw
deflection measurements (dotted line). Note: C0 is 1.4E-2 M, mercaptoethanol in water,
for this experiment.

42
 Table 2. Cantilever Parameter Estimation
for 1.4 E-2 M Mercaptoethanol in Water

Cantilever ka (M-1sec-1) kd (10-4) (sec-1) max(1015)

No. 1 0.333555 5.1484 1.27565
No. 2 0.329629 5.5826 1.32083
No. 3 0.319429 4.3686 1.35710
No. 4 0.329063 4.7588 1.26563
No. 5 0.310646 4.4155 1.21784
No. 6 0.342959 4.7652 1.16820

AVERAGE 0.32755 4.8398 1.26751

For comparison: ka is 1.3E-2 M-1 sec-1 for 1 mM CH3 (CH2)15 S/Au in ethanol [13]. G =
-RT Ln (ka/kd) = -24E3 J for this system; for 718 base pair thiolated DNA, G = -34.6E3
J [4]. max = 1.3E13 for 20-mer DNA [14].

Gauss-Markov Model Simulation

For our problem we chose to use G (t ) as an unknown but constant parameter

( G& (t ) = 0 ) and the nonlinear deflection and known temperature measurement, that is,
defining

Eq. 33 G (t ) = G (t 1) + w(t 1)

with cantilever array measurement

Eq. 34 yl (t ) = d l (t ) + (t ) + vl (t )

where G (t ) is the free energy at the surface, d l is the displacement of the lth -lever,
is the temperature of the fluid medium and w, vl are the additive, zero-mean, Gaussian
noise processes with covariances, Rww and Rvv R L L with diagonals,
v2 (l); l = 1,L , L . We assume that the measurement uncertainty is uncorrelated
producing the diagonal matrix. Each of the cantilevers has a different modulus associated
with its structure creating a set of Stoney equation parameters, l , one representing
each of the individual lever dynamics. A typical set of cantilever simulation data is
shown in Fig. 35 where we used a 20 dB SNR defined by

43
 Dtrue
2
(l )
Eq. 35 SNRl = ; l = 1,L , L
v (l )
2

where Dtrue is the true deflection available from the Gauss-Markov simulation and
Dtrue
2
is its variance. Once the noisy deflection measurements are synthesized, then the
processors are applied to extract the true deflections. We chose to evaluate two
methods: (i) smoothing processor; and (2) model-based processor. The smoothing
processor is simply a running window average that is equivalent to a low-pass filtering
operation. This is an example of a typical approach taken by scientists in this field. The
MBP is the optimal (approximately) solution to this problem. Note that we used the
third party toolbox, SSPACK_PC, a commercial model-based signal processing
package[15].

4000 4000
Deflection (nm)

2000 2000

0 0

-2000 -2000

0 5000 10000 0 5000 10000

Time (sec) Time (sec)

4000 4000
Deflection (nm)

2000 2000

0 0

-2000 -2000

0 5000 10000 0 5000 10000

Time (sec) Time (sec)

4000
4000
Deflection (nm)

2000 2000

0 0

-2000 -2000

0 5000 10000 0 5000 10000

Time (sec) Time (sec)

Fig. 35: Noisy Cantilever Deflection Measurement Gauss Markov Simulation (-20 dB
SNR).

Model-Based Processor Design

The MBP was developed using the cantilever measurement model of Eq. 34 with
the average parameter estimates of Table 2. In this section we discuss the application of
the MBP to simulated cantilever array measurement data at a variety of SNRs. We first
show the results of the MBP design for the 20 dB SNR case and then summarize the
results at a variety of SNRs to evaluate its overall performance. Note that even though the
data were simulated using the recorded temperature, it was removed (subtracted) from the
final estimates in calculating the corresponding error statistics. The MBP design
discussed is based not only on the average parameter estimates (see Table 3), but also

44
using a smoothed temperature data estimate, (t ) , to be more realistic in performance
evaluation. We expect this processor to provide an outstanding performance, once tuned
[10]. In fact, the results of applying to be 20 dB data indicating an approximate optimal
performance, since the underlying prediction errors or innovations associated with each
lever are statistically zero-mean and white as shown in Fig. 36 and Table 3. Table 3 also
shows the aggregated weighted-sum squared residual (WSSR) statistic indicating optimal
performance. The results of processing the 20 dB SNR deflection data are shown in Fig.
37 for each lever. The results are shown by comparing the true (synthesized) deflection
compared to the smoothed and MBP estimates. It is clear from the figure that the MBP
performs extremely well.
1.2 1.2
Cantilever No. 1 Cantilever No. 2

Norm. Corr.
Norm. Corr.

0.8 0.8

0.4 0.4

0 0
0 2000 4000 6000 0 2000 4000 6000

1.2 1.2
Cantilever No. 3 Cantilever No. 4
Norm. Corr.

Norm. Corr.
0.8 0.8

0.4 0.4

0 0

0 2000 4000 6000 0 2000 4000 6000

1.2 1.2
Cantilever No. 5 Cantilever No. 6
Norm. Corr.

Norm. Corr.

0.8 0.8

0.4 0.4

0 0
0 2000 4000 6000 0 2000 4000 6000

6200

6000
Rho

5800
WSSR Statistic
5600
0 2000 4000 6000 8000 10000

Fig. 36: Optimality Tests: Cantilever Array Zero-Mean/Whiteness and WSSR Statistic.

45
 Table 3. MBP Optimality Zero-Mean/Whiteness Test Results

Cantilever Mean Bound Whiteness (% out) Pass/Fail

No. 1 7.64 16.7 4.25 P
No. 2 7.81 16.6 4.48 P
No. 3 11.6 16.9 4.39 P
No. 4 8.98 16.5 3.59 P
No. 5 1.93 16.8 4.51 P
No. 6 4.16 16.6 4.34 P

WSSR < Threshold P

800
Cantilever No. 1 Cantilever No. 2
1000
600 Smoother
Smoother
Deflection

Deflection

400
True
500 True
200

MBP MBP
0 0

0 2000 4000 6000 8000 0 2000 4000 6000 8000

Time (sec) Time (sec)

1000 800
Cantilever No. 3 Cantilever No. 4
800
600
Smoother
600
Smoother
Deflection

Deflection

400
True True
400
200
200 MBP
MBP 0
0

0 2000 4000 6000 8000 0 2000 4000 6000 8000

Time (sec) Time (sec)
1000
1000
Cantilever No. 5 Cantilever No. 6
800 800

600
600 Smoother Smoother
Deflection

Deflection

400 True
400 True
200 200
MBP
0 MBP 0

0 2000 4000 6000 8000 0 2000 4000 6000 8000

Time (sec) Time (sec)

Fig. 37: MBP and Smoother Enhancement of the Noisy (-20 dB SNR) Synthesized
Deflection Measurements: True Deflection (with temperature data), Smoothed and MBP
Estimates.

46
 Next we investigate the overall performance of both the standard
smoother/averager and the MBP on synthesized data sets. In order to quantitatively
evaluate the performance of the processors individually we calculate the residual
deflection errors defined by

Eq. 36 d%l (t ) := dtruel (t ) dl (t ) ,

where dtruel is the true (noise free) deflection at the lth -lever and dl is the filtered or
estimated deflections as shown for the 20 dB case in Fig. 37. Once the error is estimated
at each lever then its associated statistics mean and variance can be calculated and used
for further analysis. We define the processor output gain as the final metric given by

d2 (l)
Eq. 37 SNRout (l) := true
; l =1,L ,L .
d2% (l)

This ratio represents the enhancement provided by each processor. The smaller the
residual error variance ( d2% ), the higher the SNRout , providing a reasonable metric.
Averaging these statistics over the cantilever array gives a feel for overall processing
gain.

We performed a set of simulations at 0, -20, -40 dB SNRin with the typical run
outputs shown in the previous figures for the 20 dB case. The results for each run are
summarized in Table 4. Here we see that output SNR produced by both processors are
respectable with a significant gain in enhancement; however, it is clear that the MBP
consistently demonstrates superior performance with an overall average enhancement of
80 dB and enhancement gain over the Smoother of 38, 43 and 60 dB, respectively. The
MBP is insensitive at these SNRs to the measurement noise variance changes yielding
identical performance at each level. By contrast, the Smoother performance clearly
deteriorates as the input SNR decreases. The MBP performance will also deteriorate with
decreases in input SNR, but not at realistic experimental levels. This demonstrates
outstanding performance for multi-channel cantilever arrays on these simulated data sets.
Next we apply the process to the measured array data.

47
 Table 4. MBP/Smoother Performance Analysis

SNRin(dB) 0 0 -20 -20 -40 -40

Cantilever SNRo(Smt) SNRo(MBP) SNRo(Smt) SNRo(MBP) SNRo(Smt) SNRo(MBP)

No. 1 44.0083 79.9439 36.5903 79.9439 17.6070 79.9440

No. 2 46.0711 87.2620 40.2018 87.2620 26.3940 87.2620
No. 3 44.9283 83.3775 44.5746 83.3775 23.5760 83.3770
No. 4 45.2004 80.3783 37.9417 80.3783 17.3240 80.3780
No. 5 44.8704 83.8749 39.8227 83.8749 22.9630 83.8750
No. 6 45.9690 82.4926 40.0839 82.4926 24.0330 82.4930

AVG 45.1746 82.8882 39.8691 82.8882 21.983 82.888

AvgGain +37.7136 + 43.0191 +60.905

5.4 Model-based processor application to multichannel cantilever data

In this section we developed MBP for two cases: (1) average deflection data and
an averaged cantilever model; and (2) multi-channel deflection data. First we averaged
the 6 cantilever deflection data to obtain the average deflection response over the entire
array and then designed the MBP obtaining the parameter estimates first and then
applying the processor with those parameters. Next we developed the multi-channel
processor using the same approach by estimating the individual cantilever parameters,
taking their average and calculating the appropriate moduli for the multi-channel model.
In both cases we performed simulations first to develop the optimal MBP.

The average cantilever model was developed by performing the model-based

parameter estimation (see Fig. 9c) obtaining the adsorption, desorption and free energy
and then developing the corresponding Gauss-Markov simulation model and
corresponding MBP. The raw and simulated data (0 dB) are shown in Fig. 38a and 38b. It
is clear that the measurement noise severely distorts the desired deflection signal. The
enhanced MBP output (deflection) enhancement is shown in Fig. 38c along with the
corresponding optimality tests in 38d where we see that the performance of the processor
is indeed optimal since the corresponding innovations are zero-mean and white. Next we
applied the MBP to the actual deflection and temperature profile data shown in Fig. 39a
and 39b after tuning the noise covariance parameters ( Rvv ) with the results shown in
Fig. 39d. Here we see that the MBP is capable of tracking the averaged cantilever
deflection data quite well; however, the performance is suboptimal, since the innovation,
although quite small, are not white. Generally, the MBP performance for this data is quite
good.

48
 3000
(a) TRUE MEASUREMENT (b) RAW MEASUREMENT (snr=0dB)

1000
1000
Amplitude

Amplitude
0
500

(Pct = 5.2

0 -2000
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
Time Time

1500 1
(c) (d) OPTIMALITY
MBP OUTPUT
ZERO-MEAN/WHITENESS TEST
ZMEAN: 7.13<8.78
1000 0.5

Norm. Corr.
Amplitude

WTEST: Pct = 2.30

500 0

0 -0.5
0 500 1000 1500 2000 2500 0 500 1000 1500
Time
Lag Time

Fig. 38: MBP Design for Simulated Average Cantilever Data (0 dB SNR): (a) True
measurement. (b) Simulated Gauss-Markov deflection measurement with bounds. (c)
Enhanced deflection (MBP output). (d) Optimality tests (zero-mean/whiteness) results.

1500 3

(a) RAW DEFLECTION DATA (b) RAW TEMPERATURE DATA

1000 2
Temperature (o)
Deflection (nm)

500 1

0 0

-500 -1
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Time (sec) Time (sec)

1200
PARAMETER ESTIMATION MODEL-BASED ENHANCEMENT
1200

(c) (d)
Model-Based
Estimate
800
Deflection (nm)

800
Model-Based
Deflection (nm)

Estimate
Raw Cantilever Data
Raw Cantilever Data 400
400
Innovation

Parameter Estimates: 0
Ka = 0.38300
0 Kd = 0.00039
Gmax = 1.07100e+015
-200
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Time (sec) Time (sec)

Fig. 39: Average Cantilever MBP Application to Measured Deflection Data. (a) Raw
deflection data. (b) Raw temperature profile data. (c) Model-based parameter estimation
fit and parameters. (d) Model-based enhancement including temperature.

49
 Next we developed the MBP for the multi-channel case using the same model-
based approach: simulation and application. We used the average model parameters
developed over the entire cantilever array data set using the nonlinear least-squares
model-based parameter estimator and applied it to the raw cantilever data to investigate
its performance. We again used the MBP with the free energy as our piecewise constant
parameter (state) and the nonlinear cantilever array model with 6 elements. We used the
smoothed estimate of the temperature profile in our estimator as in the synthesized data
case. We also used the estimated moduli in the Stoney equations as well as the average
parameter estimates as before. The results are shown in Fig. 40 where we see the raw
measured cantilever data, MBP estimates and the corresponding errors or innovations.
Since the innovations are not zero-mean and white, the processor is not optimal;
however, the results are quite reasonable except for the systematic bias error (constant) in
the estimate. The dynamics appear to be captured by the model especially in cantilever 5.
The offset can be adjusted by selecting various combinations of elements in the Rvv
measurement noise covariance matrix, but this may be better suited to an adaptive
implementation of the processor which will be pursued in future work. From the figure
we note that the dynamics of the individual levers (on-set and off-sets) are quite close to
the expected. This concludes the application to the data.

Cantilever No. 1 Cantilever No. 2

1000 1000
Deflection (nm)

Data
Deflection (nm)

Data
500 Innovation 500 MBP
MBP Innovation
0
0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Time (sec) Time (sec)

Cantilever No. 3 Cantilever No. 4

1000 1000
Deflection (nm)

Deflection (nm)

Data
Data
500 500 Innovation
Innovation
MBP MBP

0 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Time (sec) Time (sec)

Cantilever No. 5 Cantilever No. 6

1000 1000
Deflection (nm)

Deflection (nm)

Data
Data

500 500
Innovation Innovation
MBP MBP

0 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Time (sec) Time (sec)

Fig. 40: MBP of Experimental Cantilever Array Data: Raw, Enhanced (MBP) and
Residual (Innovation) results for each lever.

50
6. Summary

We have developed a successful model-based approach to the micro-cantilever

array signal enhancement problem. Experimental data were obtained from a multi-
cantilever detection system, then quantitatively analyzed with mathematical tools from
physical chemistry and solid state physics. We incorporated the results of this analysis
directly into the generic signal processing approach. A proof-of-concept solution was
created to parameterize our theoretical model, enabling us to test an average model. This
model was then used to develop the MBP for enhancing noisy cantilever measurements.
We investigated both data averaged over the array and the multi-channel cases.

Through simulations with additive Gaussian noise at SNR of 0, -20, and 40 dB,
we demonstrated the ability of the processor to extract the cantilever deflection response
with a large improvement in signal gain (~ 80 dB). We compared the processor to that of
an averager (low-pass filter) comparing their performance. At each of these SNR the
MBP demonstrated superior performance with an overall average processing gain of ~
40-60 dB over the averager. Finally, we applied the MBP to noisy, smoothed (averaged)
cantilever data and demonstrated that the processor could perform quite well except for a
bias error, which is easily corrected.

To demonstrate the full utility of the MBP for chemical sensing of low levels of
signature chemicals associated with nuclear processing, necessary next steps are (1)
verify the physical models used in this study for a variety of solvents and target
molecules (this data has already been obtained as part of this study) (2) make use of
control levers, and (3) extend the experimental library to include low concentrations of
chemical targets of practical interest for sensing nuclear processes.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the contributions of LLNL scientists Brad Hart,
Brandon Weeks, and Steve Velsko, who provided guidance on chemical sensing and
cantilever experimental protocols, and of DHS Fellow Asher Sinensky, who performed
important measurements and analysis for this project. Funding came from the NA-22
Exploratory Research Program managed by David Dye. This work was performed under
the auspices of the Department of Energy by the University of California/Lawrence
Livermore National Laboratory under contract W-7405-ENG-48.

51
7. Appendix A: Derivation of adsorption and desorption kinetics

In this appendix, a different set of variables is used than was used in previous sections of
the report. These parameters are more traditional in the field of kinetics, and so may be
more familiar to some readers. In order to compare equations in this appendix with
equations elsewhere in the report, use the following equivalencies:

ca = max
= /max

Consider a surface having a finite number of adsorption sites immersed in a gas or

solution of potential sorbants. At equilibrium, the rate of sorption and desorption are
equal:

ka c0 (1-e) = kd ca e (A1)

where kf is the forward rate constant (adsorption), c0 is the solution concentration, e is

the fraction of total adsorption sites occupied at equilibrium for this particular
concentration, kr is the reverse rate constant (desorption), and ca is the total areal
concentration of adsorption sites.

Solving for , we obtain the form of the traditional Langmuir isotherm:

e = b c0 / (1+b c0) (A2)

where

b = ka / kd ca (A3)

Ordinarily, b is defined as b0exp(-H/RT), where H is the enthalpy of adsorption.

Using the standard Arrhenius formalism for the ks, we find

Aa e-Ea/RT/ ca Ad e-Ed/RT = b0e-H/RT (A4)

At which point it is obvious that b0 = Af /caAr and H = Ea Ed. The Langmuir isotherm
has the property that equilibrium fractional surface coverage is linear with concentration
at low concentration and reaches an asymptote of one at high concentration. Also note
that the fractional coverage as a function of concentration is a function of temperature.

Now consider the rate of change in surface coverage away from equilibrium. The net
change in fractional surface coverage is the difference between the forward and reverse
rates:

d/dt = ka c0 (1-) kd ca (A5)

52
Note that the net rate approaches zero as the surface concentration approaches its
equilibrium value, when the forward and reverse rates by definition must become equal.
This implies that one can neglect desorption when deriving ka only for very low surface
coverages. One can factor the equation to get a relationship describing the rate as a
function of the distance from equilibrium:

d/dt = kac0(1-- /bc0) = kac0(1- (1+1/bcs)) = kac0(1- (bc0+1)/bc0) (A6)

Comparing to Eq. A2, we find

d/dt = ka c0 (1-/e) (A7)

The rate goes to zero when = e. Although Eq. A7 is factored to approach equilibrium
from the adsorption side, one can also factor Eq. A5 to approach equilibrium from the
desorption side:

d/dt = kd ca (- +(1-)bc0) (A8)

If c0 = 0, one gets the trivial result that

d/dt = -kd ca (A9)

Note that all adsorption terms have fallen out automatically. Now consider the case
where the concentration has been reduced (but not eliminated) so that the surface
coverage is above its equilibrium value. The resulting approach to the new equilibrium
value is given by

d/dt = kdcabc0 (1-(1+bc0)/bc0) = kdcabc0 (1-/e) (A10)

Note that d/dt is negative in this form because /e is greater than one. One can
eliminate bc0 in favor of e by first inverting Eq. A2 to give

bc0 = e/(1-e) (A11)

and then substituting Eq. A11 into the next to last equality, giving

d/dt = kdca (1-/e) e/(1-e) = -kdca (-e)/(1-e) (A12)

One observes an effective first-order desorption rate constant equal to the zero-
concentration value divided by the factor (1-e).

Now consider the case where the apparent desorption order is not unity, which we
observed for our gold-thiol system. This could be due to a difference in absorption
energies for different sites, which would cause the activation energy and desorption rate
constants to be different (Figure A1).

53
We have shown previously [16] that a distribution of activation energies yields a
deceleratory reaction that can be fitted well with an nth-order reaction having n>1.
Formally, an nth-order reaction corresponds to a gamma distribution of frequency factors,
but distributions of frequency factors and activation energies are difficult to distinguish
over narrow temperature ranges. An important feature of the model in Figure A1 is that
the non-unity reaction order would be observed only in the desorption direction, because
the adsorption direction is characterized by a single activation energy.

Ef
Er,i

reaction coordinate
Figure A1. Conceptual model giving a distribution of reaction constants for desorption,
which is reflected in an effective reaction order greater than one.

54
REFERENCES:

1. Hagan, M.F., A. Majumdar, and A.K. Chakraborty, Nanomechanical forces

generated by surface grafted DNA. Journal of Physical Chemistry B, 2002.
106(39): p. 10163-10173.
2. Halperin, A., A. Buhot, and E.B. Zhulina, Sensitivity, specificity, and the
hybridization kinetics of DNA chips. Biophysical Journal, 2004. 86(2): p. 718-
730.
3. Lavrik, N.V., et al., Gold Nano-structures for transduction of biomolecular
interactions into micometer scale movements. Biomedical Microdevices, 2001.
3(1): p. 35-44.
4. Marie, R., et al., Adsorption kinetics and mechanical properties of thiol-modified
DNA-oligos on gold investigated by microcantilever sensors. Ultramicroscopy,
2002. 91(1-4): p. 29-36.
5. Georgiadis R., Peterlinz K.P., and Peterson A.W. Quantitative Measurements and
Modeling of Kinetics in Nucleic Acid Monolayer Films using SRP Spectroscopy.
J. Am. Chem. Soc. 2000. 122, 3166-3173.
6. Karpovich, D.S. and G.J. Blanchard, Direct measurement of the adsorption kinetics
of Alkanthiolate self-assembled monolayers on mirocrystalline gold surfaces.
Langmuir, 1994. 10(10): p. 3315-3322.
7. Roark and Young, Formulas for Stress and Strain. p. Sect. 7.2.
8. Sader, J.E., Surface stress induced deflections of cantilever plates with
applications to the atomic force microscope: Rectangular plates. Journal of
Applied Physics, 2001. 89(5): p. 2911-2921.
9. Haiss, W., Surface stress of clean and adsorbate-covered solids. Reports on
Progress in Physics, 2001. 64(5): p. 591-648.
10. Candy, J., Signal Processing: The Model-Based Approach. 1986, New York:
McGraw-Hill.
11. Trees, H.V., Detection, Estimation and Modulation Theory. Vol. 1, New York:
Wiley.
12. Gill, P., W. Murray, and M. Wright, Practical Optimization. 1981, New York:
Academic Press.
13. Peterlinz, K.A. and R. Georgiadis, In situ kinetics of self assembly by surface
plasmon resonance spectroscopy. Langmuir, 1996. 12: p. 4731-4740.
14. McKendry, R., et al., Multiple label-free biodetection and quantitative DNA-
binding assays on a nanomechanical cantilever array. Proceedings of the
National Academy of Sciences of the United States of America, 2002. 99(15): p.
9783-9788.
15. Candy, J. and P. Candy, SSPACK_PC, a model-based signal processing package
on personal computers. 1993, DSP Applications.
16. Burnham, A.K. and R.L. Braun, Global kinetic analysis of complex materials.
Energy and Fuels, 1999. 13(1): p. 1-22.

55