I

1130

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 11, NOVEMBER 1991

System Identification of Electrically Coupled Smooth Muscle Cells: The Passive Electrical Properties Ping Fu and Berj L. Bardakjian, Member, IEEE

Abstract-A system model approach based on a network model is used to investigate the passive electrical properties of coupled smooth muscle cells. This approach makes use of a gradient method of optimization to estimate the passive electrical parameters directly from the magnitude of the input impedance or voltage transfer function of the network model. The need for subjective measurements of parameters and many intermediate steps involved in the analysis using the conventional signal model approach are eliminated. The coupling resistance and capacitance are estimated with sound theoretical and mathematical grounds directly from experimental data. From the simulated results using SPICE [26], it is evident that the system model approach is accurate, flegible, and reproducible. These properties grant the new approach excellent potential for future studies of drug actions on smooth muscle cells and their associated electrical coupling. Also, sensitivities of the network model with respect to its parameters can readily be obtained. This may provide new insight into the coupling mechanisms of smooth muscle cells.

I. INTRODUCTION IFFERENT types of smooth muscle cells have been shown to exhibit passive electrical properties [ 11-[6]. Since the coordination of smooth muscle cells relies heavily on electric current [7], a quantification of the passive electrical properties will elucidate the nature of the coupling mechanisms between smooth muscle cells. Excitable cell membrane is known to demonstrate both. resistive and capacitive properties. Altering these properties provides a signal which enable a cell to communicate with its surrounding. Another mode crf signaling is carried out by electrical coupling, which is formed by a protein called connexin [8], and apposing membranes in the vicinity, to provide a low resistance pathway for the flow of information among different cells. The ultrastructure of these couplings in several types of excitable tissues has been investigated using different techniques [9]-[ 141. With the structural knowledge of smooth muscle, a cable model was proposed [ 151 to analyze the passive electrical properties of smooth muscle. Since then the passive electrical properties of excitable cells have been used in

D

Manuscript received August 25, 1989; revised May 20, 1991. This work was supported by the Natural Sciences and Engineering Council of Canada (NSERC) and the Medical Research Council of Canada (MRC). The authors are with the Department of Electrical Engineering, Institute of Biomedical Engineering, and Playfair Neuroscience Unit, University of Toronto, Ont., Canada M5S 1A4. IEEE Log Number 9103286.

various branches of medicine, such as gastroenterology [ 161, cardiology [ 171, and obstetrics [6]. However, the estimation of the passive electrical parameters has been difficult, primarily due to the complexity of the smooth muscle structure and the lack of an objective technique for parameter estimation. There are two major types of modeling technique, namely the signal model and the system model. The signal model estimates the signal parameters (membrane time constant, T,, and space constant, A) using a time domain analysis of voltage response to current stimulation. The system model, however, identifies the system parameters (membrane, myoplasmic, and junctional impedances) using a frequency domain analysis of input impedance or voltage transfer function. The system model is more advantageous for the following reasons: 1) The system parameters provide a complete description of the system, whereas the signal parameters only characterize certain propagation properties of the system. 2) System models provide higher degrees of freedom. Five variables are estimated directly, compared to two in signal model. 3) Signal parameters are easily computed from system parameters, but the reverse is not necessarily true. 4) In order to identify the system parameters using a signal model, numerous intermediate steps are required. In addition, a one-to-one relationship is not guaranteed, because five system parameters are estimated from two signal parameters. The signal model analysis applied to the smooth muscle cable is similar to that used for modeling of a single neuron [18], [19]. h and T,, are subjectively estimated from voltage transients. The membrane resistance and capacitance are in turn computed indirectly using X and T,,,. Since the cable equations were developed specifically for an infinite uniform cable, the nonuniformity introduced by the anatomical separation of cells had to be neglected [22]-[25]. The assumption is true only if 1) the junctional resistance is equal to the myoplasmic resistance, and 2) the junctional time constant, rj, is much smaller than 7,. Our biological data [23] shows that both 1) and 2) are not valid. In fact, signal discontinuity of gap junctions had been observed in biological experiments previously [2]. Also, the distance between the stimulating electrode plates was required to be much smaller than the space constant. In addition, the coupling resistance and capacitance were

0018-9294/91$01.00 @ 1991 IEEE

1

FU AND BARDAKJIAN: SYSTEM IDENTIFICATION OF SMOOTH MUSCLE CELLS

..

I,

,

,

I.

..

..

I I

I

r.

I

1131

I

I I

I

unit j

I I

I

' I

I

1_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ I

I , L

I

________ _ _ _ _ _ _ _ _ _ _ _ _ _ _J

(b) Fig. I . (a) A syncytium model: two cells are connected by a gap junction. (b) A spatially discretized version o f the syncytium model: two consecutive units are shown here. Each unit is represented by a number o f compartments. Each compartment is made up of i ) a parallel combination of a membrane resistance ( R , ) and capacitance ( C M ) ,and i i ) a myoplasmic resistance ( R A ) .The gap junction is represented by a resistance ( R J )and a capacitance ( C J ) .Each node is designated by two indexes: the first is the unit number, the second is the compartment number.

estimated cnidely using some indirect means, necessitating more assumptions and approximations. This further degraded the usefulness of the signal model approach. A system model approach based on a network model 1201-[24] was developed to investigate the passive electrical properties of smooth muscle. The network model is a discretized version of the conventional cable model. The new approach eliminates the need for subjective measurements of parameters and the use of cable equations. The estimation of the membrane resistance and capacitance no longer requires the numerous intermediate steps involved in the signal model. More importantly, the coupling impedance is estimated, for the first time, with sound theoretical and mathematical grounds directly from experimental data. Moreover, the sensitivities of the network with respect to the resistances and capacitances, as well as the time constants can readily be found, which will provide new insight into the passive electrical properties of smooth muscle. 11. METHODS Derivations for the input impedance and voltage transfer function for a network model of a population of smooth muscle cells will be presented in this section. A. The Network Model The network model [ 151 consists of a number of equivalent units, which together represent a population of smooth muscle cells. Each unit is composed of a number of compartments, each of which consists of a membrane resistance (RM),membrane capacitance (C,), and myo-

plasmic resistance (RA). Adjacent cells are connected by a gap junction made up of a resistance (R,) and capacitance (C,) [Fig. l(b)]. This network is a discretized version of the previously used cable model [Fig. 1(a)]. Recent publications point to the fact that a protein called connexin constitutes the gap junctions [8], and a purely resistive junctions is suggested. However, communication between smooth muscle cells is not necessarily limited to the channels formed by connections [3]. The apposing membrane in the vicinity of the channels can play a part too. Membrane capacitance [7] can be introduced at the junctions by interdigitations [ 141 of apposing membrane. In addition, the shape of the junctional potential and current [5] obtained from intracellular recordings cannot be explained by a pure resistance. The spatial discretization of the cable model is done by dividing the cable into many equal segments. The membrane and myoplasmic resistances and the membrane capacitance are then determined depending on the number of segments used. This discretization process raises the question of how many units (M) and how many compartments ( N ) should be used. It was found that for a network model of a single neuron, which is similar to a unit in our network here, the magnitude of the input impedance and a number of electrotonic parameters were not affected significantly when more than 20 compartments were used (N > 20) [20]. Hence, similar thresholds for M and N were expected in our smooth muscle network model. SPICE [26] simulations of the network model were used to investigate the effects of discretizing the smooth muscle model with different number of M ' s and N s on the magnitude of the input impedance of the network.

I

1132

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. I I , NOVEMBER 1991

B. Derivation of the Input Impedance A spatial transform had been used ,201, [2 1, to express the input impedance of the network model of a single neuron. A similar method can be applied to the smooth muscle model. Consider a compartmental RC network in Fig. l(b) where j and k are the indexes for a unit and a node within a unit, respectively. The voltage and current relations between nodesj, k + 1 can be written as

Taking the voltage drop across the gap junction into consideration, the voltage and current expressions at node [ j + 1, 01 can be written as

I[;

-

+ 1 , 01 -I

where

61

=

1

=

RA

CI =

+ RA(GM + s C M )

(24 (2b)

GM

+ SCM

(9) where

(2c)

d, = 1

Gj

(24

G M = l/RM.

W,.k+I = HW,.,. (3) Taking the z transform of (3) with respect to the index k gives

1/Rj

YJ = GJ

(2e)

For compactness in subsequent operations, (1) is rewritten as

=

(loa)

+ JWCJ

( 1 Ob)

U = l / Y J - RA.

(10c)

Similar to steps (3) to (7), a second z transform with respect to the indexj yields V[M

[IlM

+ 1 , 01

+

1 , 01

1

(4) where

where and W,(z) is the z transform of the equally spaced sequence W,.k.Rearranging (4) gives

p2

w,(z)

=

z(zZ - H)-I

?,(I

(6)

where Z is a 2 x 2 identity matrix. Taking the inverse transform, then

z

=

[(Pz +

=

0 . 5 ( ~ 2+

- (P2 -

92)‘

42)‘I

(12a) ( 12b)

62)

~2 =

0.5[(a2 - d?)?

a2 =

xN+l

+ (Ucl

+ 4b2~21”’

-

( 12c) (124

1)xN

arbitrary voltage at the open end. Hence,

where

- d2aM

Z;,,[M + 1 , 01 =

aM+l

c2 aM PI

=

0.5(al

+ dl)

41 = 0.5[(al - dl)’

(8b)

+ 4b,~I]’”

and N = number of compartments in a unit.

Subtracting the extra gap junction from (12) gives

(8c)

z;,,w,NI

=

(aM+I

- dZaM)

c2 a,,,

_ -1 YJ

(14) ’

~

FU AND BARDAKJIAN: SYSTEM IDENTIFICATION OF SMOOTH MUSCLE CELLS

I133

C. Derivation of the Voltage Transfer Function

As mentioned earlier, in cases where the stimulating current that flows in the tissue cannot be determined from biological experiments (for example, in a Tomita’s bath [ 15]), choosing the voltage transfer function as the characterizing function in the system model is more advantageous. Our task here is to derive the explicit expressions for the approximating functions used in the optimization strategy. Unlike the input impedance function, a voltage transfer function is dimensionless. A ratio of two voltages in our network is likely to reduce the number of independent variables. The resistances and capacitances which are used as independent variables in the input impedance can no longer be used here. The rationale is clear if we consider the voltage ratio across parallel combination of resistance and capacitance. Doubling the resistance and halving the capacitance will leave the voltage ratio unchanged. Hence, the product of resistance and capacitance should be used as a parameter. This product is the time constant, one of the widely used parameters in the signal model. A new parameter set is defined as r#I = { A A , A,, 7,,,, 7 , ) . where (154

= XJ =

RJ/RM

(1 5b)

Tm =

RMCM

( 15c)

= RjCj.

(15d)

7,

For compactness of the expression, the following variables are defined: S, = XA(1

+ s7,n)

S2 = A J / ( l f S, = AJ(1

$71)

+ STm)/(l + S ~ J =) SIS2/XA

Strategy The Fletcher’s method of optimization [ 2 7 ] , [ 2 8 ] , a gradient approach, is employed for the estimation of the network parameters from smooth muscle recordings. As the name implies, in a gradient method, the gradient of the objective function with respect to each optimizing parameter determines the direction of the next move, which eventually minimizes the objective function. The major component of an objective function is a numerical measure of the difference between an approximating function and a specified function. The magnitude of the input impedance function has been used as the characterizing function of the system model of a single neuron [20]. However, since current stimulation of smooth muscle cells in a Tomita’s bath [ 151 requires two extemal stimulating electrodes, the exact amount of the injected current cannot be determined experimentally. Hence, an altemate function, the voltage transfer function, can be used as the characterizing function. 2 transforms [29] are proven to be effective and useful in expressing explicitly the approximating function based on the chosen network model. The ability to formulate explicitly the input impedance and voltage transfer function is a prerequisite for the use of any gradient method of optimization. An objective function based on the least square error can be set up as follows:

( 16a)

( 16b)

(16~)

where s is the Laplace transform variable. Based on (1 l), the voltage transfer function can be written as (stimulating at node [ I , 01):

where

D. Parameter Estimation Using an Optimization

FA = IZ,,,(~,W ( or ITF(M,N)1

E,,,

= max

(19b)

[FA(+, n) - F&)l.

( 19c)

I1

d = {RA, R M , RJ, CM, C J } and { A A , A,, 7,,,, 1 for the input impedance and voltage transfer function, respectively, F A ( d , n) and Fs(n) are the equally spaced sequences of the approximating and specifying functions, respectively, n is the sampling index, and P is the number of sample points. The approximating functions are computed from (14) and (17) for the input impedance and voltage transfer function respectively. The specified function is acquired by taking FFT’s of the recorded voltage and current transients from biological experiments. An appropriate choice of parameter space and a proper scaling of the parameters are necessary for performing the optimization processes. Linear o r nonlinear scaling can be used to ensure that the scaled model parameters are in the same order of magnitude. The model parameters should also be constrained to be positive. The gradients of the respective approximating function with respect to the scaled constrained parameters can readily be obtained by the chain rule [ 2 0 ] , [21]. 7J

- 1

1134

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 11, NOVEMBER 1991

The solution obtained is a local optimum and hence it may not be unique in a global sense. However, it is unique in a local sense, that is, in the feasible region with the constraints that the parameters are positive and that they are confined to a reasonable range of values ( f 2 order of magnitudes). Therefore, a set of model parameters estimated by the optimization method should be perturbed and then the objective function, U ( 4 , w), should be minimized again to ensure that there is no other local minimum in the vicinity of the optimal solution. Also, several different starting points should be used in the optimization process. 111. RESULTS A. Efect of “ M ” and “ N ” The values of the network components used in the simulations are listed in Tables I-A and I-B. These values are based on the effective resistances and capacitances, looking from one end of the network, being constant. Fig. 2(a) shows that when M is fixed at 3, the magnitude of the input impedance is not affected when N > 10 is used. Similarly, from Fig. 2(b), when N is fixed at 4, the magnitudes of the input impedance superimpose exactly on each other when M > 12 is used. Thus, M = 12 and N = 10 are sufficient to represent the system model of the passive electrical properties of the smooth muscle cells in this example. However, this example provides only a clue for the minimum M and N required. In practice, for each population of smooth muscle cells, M and/or N should be increased until the system model parameters stabilize. Preliminary data from canine colonic smooth muscle [22], [23] reveals that M = 12 and N = 10 are indeed sufficient in most cases.

TABLE I VALUES OF NETWORKPARAMETERS Usb.~)I N INVESTIGATINGT H ~THRESHOLD : FOR M A N D N USING SPICE. (a) M = 3 A N D N I S VARIED. (b) N = 4 A N D M I S VARIED. Parameters

N = 2

N = 4

N = IO

N = IS

R, (Q) RM (KQ) RJ (KQ) CM( P F )

660

330 IO0 3.3 0. I

I32 250 3.3

88 37 s 3.3 0.027 0. I

CJ ( P F )

so 3.3 0.2 0. I

0.1



0.04 0.1

(A)

Parameters

M=2

M = 12

R, (0)

495

82.5 400

RM (KQ) Rj (KO) CM ( P F ) CJ ( P F )

61 6.6 0.15 0.05

0.6 0.02s 0.55

M = 20

49.5

667 0.35 0.015 0.95

(B)

B. Ver$cation of the Derived Expressions for the Characterizing Functions The derived expressions for the input impedance and voltage transfer function are verified by SPICE simulation. Fig. 3(a) shows the input impedance obtained from (14) superimposing on that obtained from SPICE. The values of the parameters are those shown as “original values” in Table I11 and IV. Fig. 3(b) shows the superimposed voltage transfer function obtained from (17) on that obtained from SPICE simulation. The parameter values are listed in Table I-B (M = 12, N = 4). C. Parameter Estimation Numerical examples from SPICE simulations will be used to validate the applicability of the system approach in modeling a population of smooth muscle cells coupled by gap junctions. The sensitivities of the input impedance and voltage transfer function with respect to their parameters are also included. The values of the resistances and capacitances used in the SPICE simulations are taken from the literature on smooth muscle [ 151. Sinusoidal stimulations of various frequencies are injected into the simulated network. The resulting input impedance or voltage

Fig. 2 . The effect of using an increasing number of M and N on the magnitude of the input impedance is investigated. It is found that IZ,” I is not significantly atlected when M > 12 and N > IO. (a) M = 3, and N = 2. 4. 10. 15. N = I O is found to be Sufficient. (b) N = 4, and M = 2, 12. 20. The threshold is M = 12.

AND BARDAKJIAN: SYSTEM IDENTIFICATION OF SMOOTH MUSCLE CELLS

I I35

TABLE I1 THESYMBOLS A N D THEIR ASSOCIATED UNITSUSEDI N T H E MODELING OF A POPULATION OF SMOOTH MUSCLE CELLS. SYMBOLS WITH SMALL LETTER SUBSCRIPT REPRESENT NORMALIZED Q U A N T I T IWITH E S RESPECT TO PHYSICAL DIMENSIONS, WHILE T H O S WITH ~ . CAPITAL LETTERREPRESENT ACTUAL QUANTITIES. UNIFORM M E M B R A RESISTANCES NE ANI) CAPACITANCES A R E ASSUMED List of Symbols

Units

I = physical length of a cable I, = physical length of a gap junction A, A,

2 1

11

21

11

-

51

41

01

71

61

91

(.I1

(a)

cm cm cm’ an’ cm2 Q cm Cl cm S2 cm2 F/cni2 F/cm ms

cross sectional area of cable surface area of a compartment A, = area of a gap junction R,, = myoplasmic resistivity RI = junctional resistivity R,,, = membrane specific resistance C,,, = membrane specific capacity C, = junctional capacity r,,, = membrane time constant = R,,, . C,,, = R , . C , ins rl = junctional time constant = RI . C, = R, . CJ X = space constant cm Z , , ( s ) = input impedance n R,,, = input resistance n M = number of units N = number of compartments in a unit R/Compartment R.4 = R , , l / ( A

System identification of electrically coupled smooth muscle cells: the passive electrical properties.

A system model approach based on a network model is used to investigate the passive electrical properties of coupled smooth muscle cells. This approac...
971KB Sizes 0 Downloads 0 Views