Modelling the vascular response to sympathetic postganglionic nerve activity

This paper explores the influence of burst properties of the sympathetic nervous system on arterial contractility. Specifically, a mathematical model is constructed of the pathway from action potential generation in a sympathetic postganglionic neurone to contraction of an arterial smooth muscle cell. The differential equation model is a synthesis of models of the individual physiological processes, and is shown to be consistent with physiological data. The model is found to be unresponsive to tonic (regular) stimulation at typical frequencies recorded in sympathetic efferents. However, when stimulated at the same average frequency, but with repetitive respiratory-modulated burst patterns, it produces marked contractions. Moreover, the contractile force produced is found to be highly dependent on the number of spikes in each burst. In particular, when the model is driven by preganglionic spike trains recorded from wild-type and spontaneously hypertensive rats (which have increased spiking during each burst) the contractile force was found to be 10-fold greater in the hypertensive case. An explanation is provided in terms of the summative increased release of noradrenaline. Furthermore, the results suggest the marked effect that hypertensive spike trains had on smooth muscle cell tone can provide a significant contribution to the pathology of hypertension.


Introduction
The contractile behaviour of smooth muscle cells (SMCs) in resistance arteries is known to be controlled by the sympathetic nervous system, whose activity exhibits bursting rhythms. In particular, the activity of sympathetic postganglionic neurones innervating such cells comprises intermittent bursts of varying amplitude and frequency (Malpas, 1998). The bursting discharge drives arterial smooth muscle cell contractions, causing vasoconstriction and increases in arterial blood pressure. This activity, in both pre-and postganglionic neurones, is known to exhibit bursting that is entrained to the respiratory rhythm, due to central coupling with respiratory pattern generators (Habler et al., 1994). A quantification of how this bursting activity influences arterial smooth muscle cell contractility is important as changes to sympathetic bursting are characteristic of cardiovascular diseases, such as hypertension and heart failure. This sympathetic nerve activity is known to be elevated in the spontaneously hypertensive (SH) rat and to exhibit amplified sympathetic-respiratory coupling (Simms et al., 2009). Recently, intracellular studies in the SH rat have shown that preganglionic neurones exhibit amplified respiratory modulation in their discharge pattern (vs normotensive Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/yjtbi controls; Wistar-Kyoto rats, WKY), characterised by 2-3 more action potentials per respiratory burst . However, this amounts to a relatively small (E1 Hz) increase in average firing frequency, which is still below the threshold for a mechanical response in many sympathetic effectors (Nilsson et al., 1985;Lacroix et al., 1988;Pernow et al., 1989). Because of this, the contribution of the amplified respiratory-sympathetic coupling to vasoconstriction and the pathology of hypertension is unclear. As we shall see in this paper though, it is not the average firing frequency alone that accounts for change in SMC contractility.
Specifically, to investigate and quantify the influence of sympathetic bursts on the degree of smooth muscle cell tone, we have built a mathematical model of sympathetic transmission to a single arterial smooth muscle cell (SMC). Several sets of differential equations are coupled to provide a pathway from action potential generation in a sympathetic postganglionic neurone  to noradrenaline release from a postganglionic terminal (Yamada and Zucker, 1992), intra-SMC signalling (Li and Rinzel, 1994;Fink et al., 1999;Lemon et al., 2003;Bennett et al., 2005) and the subsequent crossbridge formation and contractile force produced by the SMC (Hai and Murphy, 1988). In particular, the latter stages of the model adopts the differential equation model described by Bennett and co-workers (Lemon et al., 2003;Bennett et al., 2005) for the behaviour of SMCs surrounding the rat tail artery. Our model begins with the neural signalling and the consequential release of noradrenaline, that drives this SMC model. We have also revised the assumptions and physiological parameters from Bennett et al. (2005) in the light of physiological data (Julien et al., 2001). This has enabled us to investigate the influence of the pattern of sympathetic postganglionic neurone activity on the contraction of arterial SMCs.
The rest of this paper is outlined as follows. Section 2 provides an overview of the model we have developed, highlighting the key assumptions made and giving reference to the appropriate literature where the individual components of the model were developed. Section 3 presents the simulation methodology and any adaptations to the models. Section 4 contains the simulation results, both those that were used to check the modelling parameters and assumptions against previously published data and also the simulation of the entire pathway. Care is taken to quantify the contractile force observed as a function of various input burst parameters. Finally, Section 5 discusses the results of our modelling in context of the literature and proposes novel hypotheses to be tested experimentally. Fig. 1 summarises the neurovascular pathway modelled-the transmission of information from a sympathetic postganglionic neurone (SPGN) to an arterial SMC. The model considers separately the mechanisms of action potential generation in a SPGN (Fig. 1A), its transmission to the distal end of the axon, causing release of noradrenaline (NA) from the axon terminal into the neuromuscular junction (Fig. 1B), the subsequent activation of the SMC triggering contraction (Fig. 1C).

Model overview
This model is novel in that it investigates the of influence sympathetic output on SMC contraction at a cellular level, and can relate sympathetic discharge patterns to contractile response. The system includes a model of a sympathetic neurone  coupled to a model of neurotransmitter release (Destexhe et al., 1994). The released noradrenaline is used to drive a model of SMC contraction (Bennett et al., 2005). By modelling the vasoconstrictor processes at this level of detail, the influence of sympathetic patterning on arterial smooth muscle cell contractility can for the first time be fully quantified and add to our understanding of increased vascular tone in pathological conditions such as hypertension.

Postganglionic action potential generation and propagation
Action potentials are generated in the sympathetic postganglionic neurone. We use a Hodgkin-Huxley style partial differential equation model of a sympathetic neurone that was reported by Briant et al. (2014) and implemented in NEURON (Carnevale and Hines, 2006). The code for this model has been deposited on ModelDB (senselab. med.yale.edu/modeldb; accession number: 151482). The model is represented schematically in Fig. 1A, and described in brief as follows.
A single un-branched axon emerges from the soma (length 500 μm and diameter of 0:5 μm, 20 segments), ending at the site of release of NA. The ion channel currents included in the model are depicted graphically in Fig. 1A. These currents are all present in the soma, with the exception of the leak (I pas ) that is present throughout the cell membrane. The axon has the Hodgkin-Huxley conductances I Na , I DR required for spike generation and also includes N-and L-type calcium currents I N and I L , because these channels are known to drive synaptic transmission at sympathetic nerve-endings (Rittenhouse and Zigmond, 1999). Spiking in the neurone model was driven by the injection of current pulses (2 nA Â 2 ms) into the soma.
The calcium current at the distal end of the axon, generated in response to action potentials arriving from the soma of the neurone model, is used to drive release of NA at the neuromuscular junction. An action potential arriving at the axon terminal of the postganglionic model activates high-threshold voltage-sensitive calcium currents I N and I L . This calcium influx locally elevates the intracellular concentration of calcium at the terminal/synapse of the axon, ½Ca 2 þ syn . Here and in what follows, square brackets ½Á are used to describe the concentration of a chemical species, with a subscript, in this case syn , indicating a particular sub-concentration.
The kinetics of this sympathetic calcium is determined by the calcium fluxes across the membrane of the synapse (I Ca ¼ I N þ I L and the calcium pumps) and by calcium buffering: Here, F d is Faraday's constant and d is the depth of the hemispherical calcium domain. I Ca is the sum of all calcium currents across the membrane I N þ I L , extracted from the distal end of the axon in the SPGN model. The pumping of calcium across the membrane occurs at a rate k t with Michaelis constant K. Calcium is buffered with timeconstant τ r to a concentration ½Ca 2 þ 1 . The calcium concentration ½Ca 2 þ syn drives the release of NA.

Noradrenaline release kinetics
Yamada and Zucker (1992) constructed a differential equation model of transmitter release from a neurone terminal, based on a model by Parnas and Parnas (1988). These kinetics are summarised in Fig. 1B and given by the following equations: Here, calcium ions are assumed to reversibly bind to a fusion protein F. Four calcium ions bind to this protein at a rate k b , changing it to its activated state F A . The reverse process has an unbinding rate k u . Destexhe et al. (1994) simplified this system by assuming that there exists an inexhaustible pool of pre-docked vesicles V, with concentration ½V, that are ready for activation. The activation of a vesicle V by an activated fusion protein F A has forward and backward rate-constants of k 1 and k 2 , respectively. The activated vesicle, V A , is then able to fuse to the membrane of the synaptic terminal, and subsequently release its contents into the neuromuscular junction. An activated vesicle V A releases N molecules of NA into the synaptic junction at a rate k 3 .
[NA] is the Fig. 1. The mathematically modelled pathway from SPGN excitation to smooth muscle cell activation. The mathematically modelled pathway from SPGN excitation to smooth muscle cell activation. (A) The model of a SPGN from Briant et al. (2014). Action potentials propagate down the SPGN axon to the postganglionic terminal. (B) This activates I L and I N , triggering the influx of calcium into the postganglionic terminal, increasing ½Ca 2 þ syn . Four molecules of intracellular calcium bind to a fusion protein, activating it (F A ). Once activated, the fusion protein can bind to, and consequently activate, a vesicle V. The activated vesicles, V A , are assumed to be pre-docked to the synaptic membrane; once activated, it immediately releases its NA contents into the cleft. (C 1 ) Released noradrenaline activates α 1 Rs on the SMC membrane, activating a G-protein (G). This G-protein, drives the hydrolysis of PIP 2 . Hydrolysed PIP2 cleaves to form IP3, which activates an IP 3 R located on the membrane of the SR. Activation of this receptor causes an efflux of Ca 2 þ from the SR (J IP3 ), increasing ½Ca 2 þ SMC . These receptors also have inactivation and activation sites for ½Ca 2 þ SMC . Fluxes of Ca 2 þ across the SR membrane also exist due to leakage (J leak ) and calcium pumps (J pump ). (C 2 ) The intra-SMC matrix contains actin (A) and myosin (M) filaments. At rest these filaments are in a detached state, A þ M. When the myosin heads are phosphorylated by calcium (M P ), they able to attach to the actin filaments, yielding the state AM P -a cross-bridge. This cross-bridge can then conduct a 'power stroke', sliding the actin filament and producing a contractile force.
concentration of the released ligand NA. The released transmitter is removed from the synaptic cleft by diffusion, degradation and reuptake, all of which are incorporated into the clearance rateconstant k h .

G-protein activation in smooth muscle cells
Released NA diffuses across the neuromuscular junction, and binds to α 1 -adrenoreceptors (α 1 R) on the SMC membrane ( Fig. 1C 1 ). This is a G q=11 -protein coupled receptor. On binding NA the α 1 R undergoes a conformation change, activating the intracellularly coupled G-protein (G q=11 ), which detaches from the receptor, dissociating into subunits. This detached/activated G-protein triggers a second messenger cascade within the SMC, leading to a contractile response. A model of the activation, phosphorylation and internalisation of a G-protein coupled receptor has previously been described by Bennett and co-workers and is used here for modelling intra-SMC signalling (Lemon et al., 2003;Bennett et al., 2005). The specific equations model the activation of an α 1 R following binding of NA, and the subsequent decoupling of a G-protein: Here, ½R S and ½R S P are the numbers of unphosphorylated and phosphorylated surface receptors, respectively. These equations assume that the binding of the ligand NA is rapid; the implications of which were formalised by Wagner and Keizer (1994). K 1 and K 2 are dissociation constants and ½R T is the total number of receptors.
The parameter ζ ¼ 0:85 takes into account the assumption that only 85% of receptors will be participating in G-protein signalling. The parameters k p , k e , k r are the rates of receptor phosphorylation, internalisation and recycling, respectively. For a full derivation of the receptor equations, see Lemon et al. (2003).
Only the unphosphorylated surface receptors ½R S are able to bind NA and activate the G-protein cascade that leads to SMC contraction. As formalised by Lemon et al. (2003), the concentration of activated Gprotein ½G satisfies an equation of the form Here, the rate-constants k a and k d are of G-protein activation and deactivation, respectively, and ½G T is the total number of G-protein molecules. The constant δ represents the background contribution of NA-unbound, unphosphorylated surface receptors to G-protein activation. The function p r is the main contribution to G-protein activity, representing the activation of G-protein molecules by NA-bound unphosphorylated surface receptors. The form of this function is given by 2.4. IP 3 -induced Ca 2 þ release The activated G-protein triggers a protein cascade, resulting in the release of Ca 2 þ from intracellular stores within the SMC. The kinetics of this process are shown in Fig. 1C 1 and can be described as follows.
An inactive G-protein G Á GDP can bind to a NA-bound α 1 R. Once NA binds to the α 1 R then guanosine diphosphate (GDP) is exchanged for guanosine-5 0 -triphosphate (GTP). The G-protein dissociates from the receptor to form the complex G Á GTP. Once dissociated, the G-protein disassembles into subunits. The active subunit is (G α Á GTP). This activated G-protein is able to bind to phospholipase C (PLC). This G-protein-PLC complex can then hydrolyse PIP 2 (that resides in the SMC membrane) to inositol 1,4,5-trisphosphate (IP 3 . This process is catalysed by Ca 2þ . We have adopted the model of this G-protein cascade derived by Bennett and co-workers (Lemon et al., 2003;Bennett et al., 2005). Here, the active G-protein complex G α Á GTP is denoted simply by G. The rate of hydrolysis of PIP 2 is r h ½PIP 2 , where [PIP 2 ] is the total number of PIP 2 molecules in the cell membrane and r h is the rate at which one molecule is hydrolysed into IP 3 , given by Here, α is an 'effective signal gain parameter', and K C is the dissociation constant for the Ca 2þ binding site on the PLC molecule. ½Ca 2 þ SMC is the intra-SMC concentration of Ca 2 þ . The rate of change of IP 3 is determined by the rate at which PIP 2 is hydrolysed and IP 3 is degraded, replenishing PIP 2 levels in the cell membrane. The means by which PIP 2 is replenished can be summarised by the degradation of IP 3 to an intermediate phospholipid, which is then phosphorylated and returned to the cell surface. These two processes can be modelled as follows: As in Lemon et al. (2003) and Bennett et al. (2005), we require IP 3 as a concentration for subsequent equations. Therefore, [IP 3 ] is the concentration of IP 3 , [PIP 2 ] is the number of PIP 2 molecules and the factor γ (Avogadro's number, N a Â cell volume, C V ) converts concentration into particle number. [ðPIP 2 Þ T ] is the total pool of phospholipid to which IP 3 is returned. The replenishment of PIP 2 levels in the membrane assumes degradation of IP 3 to the intermediate phospholipid at a rate k deg , and then phosphorylation of this phospholipid to PIP 2 at a rate r r . The IP 3 molecule activates IP 3 receptors (IP 3 R) on the membrane of the sarcoplasmic reticulum (SR) within the SMC. Activation of this receptor causes the release of calcium into the cytosol of the SMC, increasing the intra-SMC concentration of Ca 2 þ , ½Ca 2 þ SMC . A twovariable model describing the calcium fluxes across an SR membrane was first described by Li and Rinzel (1994), and employed by Fink et al. (1999). The ½Ca 2 þ SMC dynamics are governed by fluxes across the SR membrane from a leak flux J leak , a calcium pump flux J pump , and a receptor flux generated by IP 3 channels J IP3 . Released calcium is buffered at a constant rate β, giving The additional variable h is a pseudo-gating variable for the IP 3 R channels. The three fluxes are given by

Smooth muscle cell contraction
The basis of contraction of an SMC involves the binding of phosphorylated myosin to actin, forming a 'cross-bridge' (Fig. 1C 2 ). The kinetic model we use to describe this process is again that used by Bennett et al. (2005), which is a reduction of the original model by Hai and Murphy (1988).
In the model, the states of actin A and myosin M are considered. These can be phosphorylated ( p ), and bound to each other. Free Ca 2 þ molecules are able to trigger a process leading to phosphorylation of the myosin head (M P ), cocking it. This process happens at a ratê k 1 ¼k m ½Ca 2 þ 4 . Phosphorylation of M is reversible, with an unbinding ratek 2 . Once the myosin head is cocked it is able to attach to actin and generate a contractile force. The cocked head 'latches' to the actin filament at a rate constantk 3 , with a de-latching ratek 4 . The latching of a cocked myosin head M P to an actin filament forms a cross-bridge AM P . The latched and cocked cross-bridge may now perform a contractile power-stroke, using the energy stored in the myosin head as adenosine diphosphate (ADP).
By assuming that there is a total, conserved concentration of myosin, ½M T ¼ ½Mþ½M P þ½AM P , it is possible to derive the simplified kinetics Finally, the contractile force produced by the SMC is assumed to be proportional to the concentration ½AM P of cross-bridges.

Model adaptation and simulation
Here we outline the parameter values of the model, any adaptations we have made to the equations or parameters, and details of simulations.

Noradrenaline release kinetic parameters
The parameter values for NA release kinetics are given in Table 1, together with references supporting the chosen value. The parameter values for the NA release model-Eqs. (1)-(4)-are as in the original model of transmitter release by Destexhe et al. (1994). The only deviation is the value of the parameter k h , which encapsulates the rate of clearance, degradation and reuptake of the released transmitter. Carbon fibre studies of NA release onto the rat tail artery, have revealed that it is cleared by neuronal reuptake and diffusion, while extraneuronal uptake is negligible (Stjarne et al., 1994). Since diffusion of NA is known to be restricted, due to the dense plexus of connective tissue surrounding release sites, we set k h ¼ 0:003 ms À 1 (Gonon et al., 1993b). We found that with this parameter value, stimulation of the SPGN model with a single pulse, resulted in a clearance of released NA that took % 1:5 s, fitting carbon-fibre observations (Stjarne et al., 1994;Stjarne, 2000). Furthermore, stimulation of the model with action potentials at 12 Hz produced a NA concentration of 100 nM (see Fig. 4). This fit experimental data of Gonon et al. (1993a), which found the concentration to reach 120 nM.

G-protein cascade in smooth muscle cells
The parameters of Eqs. (5)-(11) are given in Tables 1 and 2. All parameter values are as in the original model of a G-protein cascade following activation of a P 2 Y 2 purinergic receptor (Lemon et al., 2003;Bennett et al., 2005). However, we have used a timescale argument to reduce Eqs. (5)-(13). This also made the system less numerically stiff. Examination of the time constants reveals that receptor and G-protein dynamics occur on a timescale that is at least one order of magnitude faster than all other dynamic features of the model. Therefore to make the model less numerically stiff, and to reduce the dimension of the model we suppose that the receptor and G-protein dynamics are instantaneous and therefore these concentrations are replaced by their steady-state values. In particular, for a given noradrenaline concentration ½NA, Eqs. (5)-(7) can be replaced by the linear system where lim t-1 We can then determine the steady state of activated G-protein via This steady-state ½G can then be used in Eqs. (9)-(11) in place of the dynamic variable ½G. Fig. 2C validates this steady-state assumption. The intracellular calcium concentration ½Ca 2 þ SMC in response to ½NA concentrations of 0:1 μM and 1 μM closely matches those of the original system by Bennett et al. (2005).

IP 3 -induced Ca 2 þ release
The parameters of Eq. (12) are given in Table 2. Their values are taken from a model of Ca 2 þ dynamics following IP 3 uncaging in A7r5 cells (a cultured rat SMC line from the thoracic aorta) (Fink et al., 1999). This is the same parameter set as used in Bennett et al. (2005).

Cross-bridge formation in smooth muscle cells
The parameter values for Eqs. (14) and (15) are given in Table 2.
The parameters values are as in Bennett et al. (2005) except for the total concentration of myosin, which was set to ½M T ¼ 2:3. This gives a correspondence between the units of ½AM P and the contractile force generated by the SMC seen experimentally by Yagi et al. (1988) (see Fig. 2B). The variable ½AM P is therefore referred to as the contractile force produced by the SMC model, with units of μN.

Simulation methodology
Simulations of our single-cell model of an SPGN were performed with the programme NEURON (version 7.3; Carnevale and Hines, 2006) using a 12:5 μs time step. The response I Ca of the SPGN model was recorded at the distal end of the axon and imported into MATLAB (version 6.1) within which the rest of the model was run.
Simulations of our model were performed on a two dual-core Opteron processors 8 GB RAM node, using the computational facilities of the Advanced Computing Research Centre, University of Bristol, UK (http://www.bris.ac.uk/acrc/).

Model validation
To test the correctness of our adaptation to Bennett et al.'s (2005) model, different components of the model were directly simulated and compared with experimental data. First, the model was directly stimulated with fixed concentrations of NA to mimic exogenous bath application of NA and other α 1 R agonists ( Fig. 2A). The SMC contractile response followed a log sigmoidal relationship with [NA]. At 1 μM (log 10 ðMÞ ¼ À6) the response begins to rise, with an ED 50 ¼ 1:47 μM (log 10 ðMÞ ¼ À5:83). At % 10 μM (log 10 ðMÞ ¼ À5) a maximal contraction is reached. Experimental data of Sjoblom-Widfeldt and Nilsson (1989) is also shown. In these experiments, the active tension of mesenteric arteries in male Wistar rats in response to NA application, in the presence of extracellular calcium concentration 1.0 mM, was measured. The ED 50 of these experimental data was 1:48 μM, providing a good fit with our simulation data. In particular the model provides an excellent fit to the top part of the experimental response curves, but begins to rise at a higher NA concentration and with a steeper relationship to ½AM P . Note however that the simulations are for a single SMC, whereas the experimental data is for an arterial response (see Section 5 for further discussion).
In addition, we explored the characteristics of the actin-myosin model (Eqs. (14) and (15)) by driving it with a fixed level of ½Ca 2 þ SMC held at concentrations between 0-2 μM (Fig. 2B). Note how the simulated steady-state curve fits the experimental data of Yagi et al. (1988). Both the experimental and simulated response begin to rise at ½Ca 2 þ SMC ¼ 0:1 μM and reach a maximum at ½Ca 2 þ SMC ¼ 1 μM.
Finally, the response of ½Ca 2 þ SMC to exogenous application of NA was recorded (Fig. 2C). Fig. 2C shows the dynamic ½Ca 2 þ SMC response to ½NA ¼ 0:1 μM and 1 μM. Note how the calcium concentration experiences a delayed transient peak, followed by a steady-state plateau, as in experimental observations (Li et al., 1993). These simulations also fit those of Bennett et al. (2005), justifying our steady-state approximation of the G-protein and receptor dynamics given in Eqs. (16) and (19).
A further test to the model was to expose it periodically to a 300 ms bolus of ½NA (2 mM) at frequencies 0.01-2 Hz (Fig. 3). The steady-state mean ½AM P and amplitude of any oscillations in ½AM P are plotted against frequency in Fig. 3A. At high frequencies, the model cell contracted and was unable to relax before the next bolus of NA occurred; the magnitude of oscillations was seen to decrease quickly with increasing stimulation frequency. Such a tonic contraction occurred at stimulation frequencies above 0.4 Hz. At higher stimulation frequencies, the mean contractile force is greater, but the oscillation amplitude reduced (Fig. 3B). These results match those found experimentally in vascular SMCs isolated from the aorta of rats (Julien et al., 2001), where SMCs only produced oscillatory contractions (in response to repetitive boluses of the α 1 R agonist phenylephrine) when the stimulation frequency was r 0:3 Hz.

Tonic stimulation of model
To determine the response of the SMC to tonic firing, the SPGN model was driven with a train of current pulses (2 nA Â 2 ms) at a range of frequencies (1-25 Hz). These frequencies are consistent with those seen in muscle vasoconstrictor (MVC) sympathetic preganglionic neurones of normotensive (Wistar-Kyoto, WKY) rats (Malpas, 1998;Stalbovskiy et al., 2014;Briant et al., 2014).
Stimulating the model with current pulses at 9 Hz (Fig. 4A) and 8 Hz (Fig. 4B) caused contraction of the smooth muscle cell model. Each pulse (Fig. 4A 1 and B 1 ) triggered an action potential in the soma, which propagated down the axon of the SPGN model (Fig. 4A 2 and B 2 ). L-and N-type calcium channels activated causing an increase in ½Ca 2 þ syn (Fig. 4A 3 and B 3 ). The concentration of activated fusion protein, ½F A was found to increase with firing frequency (Fig. 4A 4 and B 4 ). At higher stimulation frequencies, the amplitude of ½F A oscillations increases. Because F A binds to vesicles V, activating them for release, it was consequently found that both the mean concentration (data point) and pulse-to-pulse amplitude (error bars) of activated vesicles ½V A increased with stimulation frequency (Fig. 4D). Similarly, driven by the nonlinearity in ½F A it was found that both the mean and the oscillatory response of ½V A to spikes in ½Ca 2 þ syn increase non-linearly with firing frequency. Hence, as activated vesicles dock and NA is released into the synaptic cleft, it is found that ½NA increases nonlinearly with firing frequency, but is proportional to the increase in ½V A (Fig. 4E). Notice how the oscillations in ½NA were found to have negligible amplitude at high frequency, being o9 μM in magnitude.
The subsequent G-protein cascade causes a rise in the postjunctional concentration of calcium, ½Ca 2 þ SMC , as shown in Fig. 4F.
The steady-state contractile force ½AM P in response to tonic SPGN stimulation is shown in Fig. 4G. This sigmoidal relationship began rising at 7 Hz and peaked at ½AM P ¼ 1:6 μN by 12 Hz. For tonic stimulation at o 7 Hz, the contractile response of the SMC was minimal (o 1%). The simulation data had an ED 50 ¼ 7:3 Hz which is a good fit to experimental data, where the ED 50 of mesenteric arteries is between 6 and 8 Hz (Nilsson et al., 1985;Sjoblom-Widfeldt and Nilsson, 1989).

Paired-pulse stimulation of the model
To understand the relationship between firing frequency and NA release, pairs of action potentials were generated in the SPGN, with varying inter-pulse intervals (IPIs) of 40-2000 ms. We were particular interested in the response to paired pulses separated by o 100 ms, as bursting in sympathetic efferents can reach instantaneous frequencies of 4 10 Hz.
The paired action potentials propagated down the axon to drive the pre-synaptic calcium concentration and therefore trigger NA release from the axon terminal (Fig. 5). The peak response of presynaptic calcium (Ca syn -normalised to the response to first pulse) to the second pulse, was not greatly influenced by the IPI (Fig. 5A 1  and B). In particular, the amplitude of the Ca syn response to the second pulse was similar for 50 ms, 100 ms and 150 ms IPIs (Fig. 5A 1 ). The peak gain in Ca syn increased non-linearly as IPI is reduced, but this relationship was relatively gentle (compared to other variables; Fig. 5B). Fusion proteins (F) were activated (F A ) by the increase in pre-synaptic calcium. The response of F A to the second pulse (normalised by the response to the first pulse) was greatly increased as the IPI was reduced from 150 ms to 50 ms (Fig. 5A 2 ). Because of this, the gain in activated fusion proteins in response to the second pulse, could reach 2.5-3.5 times that of the first pulse for small ( o 100 ms) IPIs (Fig. 5B), much greater than the paired-pulse gain of Ca syn . The gain in activation of vesicles was similar to that of F A , being linearly related.
The paired-pulse gain for NA was also relatively large ( Fig. 5A 3 ), reaching 3-4 times that of the first pulse for short IPIs (o 100 ms; Fig. 5B). This was driven by the paired-pulse gain in F A and the slow time-course of decay of NA.
The paired-pulse gain in F A was plotted as a function of the "baseline" (see Fig. 5A 1 ) level of Ca syn (Fig. 5C). The data show that the paired-pulse gain in F A is a function of the Ca syn level at the time of arrival of the second spike; when the baseline level of Ca syn approaches 1-tantamount to the IPI being reduce-the pairedpulse gain in F A can reach 4. These data are fit with a quartic polynomial. Therefore, the large paired-pulse gain in fusion protein activation and NA release for small IPI is because F A is a quartic function of the basal level of Ca syn .

Spike-train driven model
The model was driven with experimentally recorded spike trains and for comparison tonic trains at the same average firing frequency (Fig. 6). Membrane potential was recorded in situ from muscle vasoconstrictor-like sympathetic preganglionic neurones in normotensive (WKY) and spontaneously hypertensive rats by  whole-cell patch clamping techniques Briant et al., 2014). The spike trains from cell recordings of a WKY and SH rats are shown in Fig. 6A 1 and B 1 , respectively. The discharge exhibited the typical respiratory modulation of firing, with an increase in firing occurring during inspiration (as marked by phrenic nerve activity (PNA); see Fig. 6A 1 and B 1 lower traces and shaded regions). Note that the average firing frequency of the spike trains was subthreshold for a SMC contractile response, being 1.6 Hz for the WKY and 4.6 Hz for the SH rat (see Fig. 4G). Note also that the SH spike train exhibits increased respiratory modulated bursting, with spiking during PNA that peaks at 12 Hz (vs 5 Hz in the WKY cell). This is a characteristic of sympathetic activity in hypertension (Simms et al., 2009;Briant et al., 2014). The raw voltage-traces were converted into spike trains to drive the SPGN model-at each action potential a pulse (2 ms Â 2 nA) was played into the soma, to trigger an action potential in the model. The response of the SMC model to these inputs was compared to the response to tonic inputs of the same average firing frequency. Each action potential in the spike train triggered a transient calcium current in the distal end of the axon, driving fusion protein activation and NA release the neuromuscular junction (Fig. 6A 2 and B 2 ). The amount of NA released was found to increase phasically with PNA, due to bursts of action potentials occurring during inspiration. The concentration of NA released was greater than that produced by tonic stimulation at the same average firing frequency. The SH spike train triggered the greatest release of NA (note the difference in axis) which was entrained to PNA.
The SPGN model was driven by 10 spike trains recorded from WKY (n ¼ 5) and SH (n ¼ 5) rats in Briant et al. (2014) (Fig. 6C). Each spike train consisted of 5-10 respiratory cycles of spiking data. Spiking increased during PNA (inspiration). The WKY spike trains had low (1)-(5). The variables of these equations were normalised by the single-pulse response, to measure the gain in the second pulse as a function of IPI. (A) The response of pre-synaptic calcium (Ca syn , A 1 ), activated fusion protein (F A , A 2 ) and noradrenaline (NA, A 3 ) to paired-pulse stimulation at 50 (green), 100 (blue) and 150 ms (red). The gain in the second pulse in Ca syn was similar across the IPIs. The gain in F A increased as IPI decreased. (B) Paired-pulse peak gain in Ca syn (circle), F A (diamond) and NA (square) as a function of IPI. Gain was measured as the maximal height reached in response to the second response, as a proportion of the height of the first response (gain ¼ 1). Note that at high-frequency (IPI o 100 ms) the gain in Ca syn is small (1.3-1.4) and a shallow function of IPI. Contrarily, the gain in F A is large (2.5-3.5) and a steep function of IPI. This drives greater release of NA, which exhibits an even greater gain (3-4) due to slow reuptake. (C) The relationship between the gain in activation of the fusion protein and the baseline pre-synaptic calcium level at the time of arrival of the second action potential is quartic. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) average firing frequencies (o 2:5 Hz) and were characterised by bursts of 3-5 spikes during the inspiration period. The contractile force produced by the WKY spike trains was small, being similar to the force produced by tonic stimulation at the same average firing frequency (E0.1% of the maximal force available). The SH spike trains had larger average firing frequencies (2.3-4.6 Hz), and had characteristic respiratory bursts of larger amplitude (more spikes per burst). The contractile response to the SH spike trains was 0.2-3.8% of the force available. These larger contractions were also greater than the contractions produce by tonic stimulation at the same average firing frequency. Two SH spike trains produced small contractions (E0.2% of maximum), that were still greater than tonic and WKY stimulation. These data show that SH spike trains can evoke an order-ofmagnitude increase in the contractile force produced, despite having average firing frequencies that are subthreshold for a contractile response. This is because SH spike trains exhibit amplified respiratory bursts. We therefore sought to quantify the response of the model to bursting inputs.

Bursting stimulation of model
To quantify the behaviour of SMCs to bursts, the model was stimulated with bursts of action potentials with varying properties; burst amplitude (spikes/burst, n), inter-burst frequency (f), and burst duration. These burst parameters were varied over physiological ranges to determine their effect on the contractile force ½AM P . The results are shown in Figs. 7 and 8. Fig. 7 shows the contractile force ½AM P generated quantified as a function of the inter-burst frequency f and burst amplitude n. A series of bursts characterised by fixed f and n was generated by stimulating the SPGN model with bursts of n brief 2 ms pulses in current 2 nA amplitude, incoming at a frequency f (Fig. 7A). These pulses were equally spaced across the duration of the burst (250 ms), and induced a burst of n action potentials in the SPGN model. The number of spikes n was chosen to vary between 1 and 15, covering intra-burst frequencies seen in situ for the respiratory components of bursting in normotensive (control) Wistar-Kyoto rats (n ¼2-5)  and spontaneously hypertensive (SH) rats (n¼ 6-10) . In this experimental protocol, the respiratory frequency of the rat is between 0.7 and 1 Hz. The fixed burst duration of 250 ms was chosen to mimic; first, the duration of oscillatory inputs used in studies of sympathetic transmission to the vasculature in conscious rat (Stauss and Kregel, 1996) and; second, the duration of respiratory modulated bursts exhibited by sympathetic preganglionic neurones (SPN) in situ Briant et al., 2014). The inter-burst frequencies f ranged between 0.05 and 4 Hz. This pattern of bursting was played in for ( 4 100 s), at which steadystate behaviour was reached in the model. The mean contractile force ½AM P across one period and the amplitude of any oscillations in ½AM P was recorded, and plot against (n, f) in a 3D plot (Fig. 7B 1  and B 2 ).
Note how the mean ½AM P response shows an approximately sigmoidal increase with both increasing frequency f and number of spikes/burst n (Fig. 7B 1 ). As the number of spikes per burst decreases, response against f becomes gentler and less sigmoidal. Similarly, as the frequency f at which the bursts arrives decreases, function of n becomes gentler. Fig. 7C 1 shows a heat-map for this simulation data.   were used to drive contraction in the SMC model. The WKY cell fired at an average frequency of 1.6 Hz and the SH cell at 4.6 Hz. Both recordings exhibit respiratory modulation, with increases in firing frequency entrained to PNA (lower trace; shaded region). The response of the model to these patterns was compared to tonic stimulation at the same average firing frequency. The NA released in response to the bursting spike train, was greater than tonic, for both the WKY (A 2 ) and SH (B 2 ) recording. The release of NA was entrained to respiration. The contractile response was greater in response to the bursting spike trains, but was greatest in response to the SH spike train.
(C) Steady-state contractile response of SMC model (% of maximum) to n ¼ 5 WKY spike trains (square) and n ¼ 5 SH spike trains (circle). The response to tonic stimulation at the same average firing frequency is also shown (dashed line). The colour map indicates the level of ½AM P for each (n, f). The heatmap shows that the contractile force ½AM P has a much steeper relationship with n than it does with the frequency f at which the bursts arrive. An oscillatory response of ½AM P to bursting was not seen for all values of f and n (Fig. 7B 2 ). Only bursts incoming at low frequencies (f r 0:3 Hz) and with high amplitude (n Z 9) were able to induce significant (Z 5%) oscillations in ½AM P . When bursts consisting of n 4 10 spikes where incoming at f ¼ 0:05 Hz, the SMC was able to produce ½AM P large oscillations (35% of maximal ½AM P ). For interburst frequencies f r 0:3 Hz oscillations would begin at some critical number of spikes n, and rise to a plateau, at which further addition of spikes could not achieve greater oscillatory amplitudes. In Fig. 7C 2 a heat-map is shown for the oscillatory ½AM P data. The colour of each (n, f) square indicates the peak-to-trough change in ½AM P over one period, according to the colour map.
When bursts are played in at inter-burst frequencies matching respiratory frequencies of 0.5 Hz, as seen in the rat in situ Briant et al., 2014), burst amplitude (n) greatly increases the contractile force produced at a critical value of n ¼ 8. This critical value is a burst amplitude rarely exhibited by SPN in normotensive control rats, but characteristic of respiratory bursts in SH rats Briant et al., 2014). Fig. 8 shows the contractile force ½AM P as a function of burst duration. For these data, the inter-burst frequency was fixed at f¼ 0.5 Hz which is consistent with the frequency of respiratory modulation of SNA in situ. Bursts with fixed amplitude n ¼4 and n ¼8 spikes, yielding average firing frequencies of 2 Hz and 4 Hz, respectively, were played into the SPGN model. These burst amplitudes covered average firing frequencies and burst amplitudes seen in situ in normotensive (control) WKY rats and SH rats Briant et al., 2014). As the duration of the burst was decreased, the mean contractile force increased. The contractile force is therefore dependent on the intra-burst frequency.

Discussion
The aim of this study has been to investigate the sensitivity of arterial smooth muscle cell tone to properties of the burst signals in sympathetic postganglionic neurones. To achieve this we have synthesised and adapted previously published mathematical models (Hai and Murphy, 1988;Destexhe et al., 1994;Li and Rinzel, 1994;Lemon et al., 2003;Bennett et al., 2005;Briant et al., 2014), to produce an in silico description of the pathway from action potential generation in a sympathetic postganglionic neurone to the contractile force produced by a SMC. We have for the first time been able to explore in a model the influence of sympathetic patterning on the contractile response of smooth muscle cells. The oscillatory response of ½AM P to the bursting stimulation, measured as the peak-to-trough change in ½AM P over one period (% of max). An oscillatory response was only seen for very high amplitude (large n) bursts, at very low inter-burst frequencies f. (C) Heat maps depicting the data for (B). Note that at frequencies mimicking respiratory modulation of sympathetic activity ( E0.4 Hz, as in Fig. 6), ½AM P is very sensitive to changes in burst amplitude n. The magnitude of the [AM p ] response, and any oscillations in [AM p ] produced, is also strongly dependent on f, particularly at lower inter-burst frequencies (f o 0:4 Hz).
We found that SMCs are unresponsive to tonic stimulation at average firing frequencies typical of those reported in single sympathetic fibres (Habler et al., 1994), which highlights the importance of bursting. We showed that the SMC was able to generate contractile forces in response to these bursts (compared to tonic stimulation at the same average firing frequency) and that the magnitude of this contractile force is strongly dependent on sympathetic burst properties. In particular, adding 1-2 spikes to a burst could cause the model to generate a contractile response. Given that this increase in burst amplitude is characteristic of sympathetic preganglionic neurones in the spontaneously hypertensive rat, we drove the model with such spike trains. The contractile response generated by the model increased markedly (compared to normotensive control trains) due to potentiated release of NA. These data provides the first mechanistic understanding of how increased respiratory modulation, recorded in the pre-hypertensive SH rat, causes higher vascular resistance and therefore contributes to the pathology of this disease (Simms et al., 2009;Moraes et al., 2014;Briant et al., 2014).

What causes the greater response to bursting?
The impact of irregular or bursting sympathetic stimulation on vasoconstriction has been studied in the rat mesentery (Nilsson et al., 1985), pig nasal mucosa (Lacroix et al., 1988), dog gracilis muscle and pig spleen (Pernow et al., 1989). All studies reported an increased response of the innervated vasculature to irregular or burst patterns (compared to tonic), despite average firing frequencies (all o 5 Hz) being subthreshold for a mechanical response. Our simulations show that the ability of the vasculature to respond to such low firing rates is achieved because release of NA is potentiated by grouping spikes. Stimulation of our model with spike trains, characterised by grouped stimuli that occur during the inspiratory period (Fig. 6), produced greater contractile forces despite the SMC model being unresponsive to such a low average firing frequency. Julien et al. (2001) stimulated SMCs with phenylephrine and attributed the frequency-dependent response to NA kinetics, rather than to post-junctional intra-SMC signal processing. The present study supports this, showing that the kinetics of NA release determines the amplified response of the SMC to bursting. In cardiac myocytes, systolic calcium concentration (the amplitude of Ca 2 þ transients) depends cubically on SR content (Trafford et al., 2001). This highly non-linear relationship means that small changes in SR content markedly influence myocyte contractility. Similarly in the present study, the relationship between firing frequency and the amplitude of pre-synaptic calcium transients is highly non-linear. As the frequency increases, the transients achieve greater amplitudes, potentiating NA release (Figs. 4 and 5). This occurs via a non-linearity with activation of fusion protein; the gain in fusion protein activation was seen to be strongly dependent on the intra-burst frequency. The pre-synaptic calcium dynamics are therefore causing the preferential response to bursting; this is supported by experimental data of sympathetic stimulation of small mesenteric arteries (Sjoblom-Widfeldt and Nilsson, 1989). In the study, a moderate reduction in extracellular calcium greatly changed the frequency-response curve, a shift that was due to pre-junctional mechanisms (such as altered NA release from nerveendings).

Sympathetic bursts are not unitary events
The simulations from the present study show that the properties of a burst are important determinants of the contractile response of SMCs (Figs. 7 and 8). Analysis of muscle sympathetic nerve activity (mSNA) in humans has largely considered bursts in mSNA as unitary events-measuring bursts/min and ignoring the temporal spacing and amplitude of these burst events (Macefield, 2013). This assumption has recently been challenged by the work of Fairfax et al. (2013). They show that individual cardiac-triggered bursts of mSNA can dynamically regulate forearm vascular conductance (FVC), and that both the incoming pattern and magnitude of the bursts are able to decrease FVC by an α-adrenoreceptor mediated mechanism. Our data not only corroborates these findings, highlighting the need to move away from considering bursts as unitary events, but quantifies in detail the influence of burst amplitude, pattern and duration on the contractile response of SMCs.
The nature of the contractile response of a single SMC to periodic phenylephrine application is known to be frequency-dependent. Julien et al. (2001) studied SMCs isolated from the aorta of fourweek old rats, and demonstrated that they tonically contract when the frequency of application was 4 0:1 Hz, but could produce periodic contractions at lower frequencies. Concordantly, our model shows that only slow (o0:4 Hz) bolus application of NA to the SMC model may result in a significant periodic constriction and dilation of a single SMC (Fig. 3B); at higher frequencies of stimulation the model exhibited a maintained contraction. It has been shown that the sympathetic nervous system (SNS) is too 'sluggish' to transmit frequencies higher than 1 Hz in rat (see review Julien et al., 2001). Stimulation of the sympathetic nerve with bursts of fixed amplitude (n¼6) and duration (250 ms), revealed that the peripheral SNS in conscious rats can only transmit oscillations to the vasculature at frequencies o1 Hz-any higher frequencies simply contribute to the mean level of vasoconstriction (Stauss and Kregel, 1996). Contrarily, Fairfax et al. (2013) reported that the absence of a mSNA burst in a cardiac cycle is The bursts had a fixed inter-burst frequency of 0.5 Hz, to mimic respiratory modulation of bursting seen in situ . Furthermore, burst amplitudes of n ¼ 4 and n ¼ 8 were chosen to match the burst amplitudes recorded in situ in normotensive (control) WKY rats and spontaneously hypertensive rats, respectively . These patterns achieved average firing rates of 2 Hz and 4 Hz, approximately what was reported in the two strains. These spikes were then equally spread over a burst of duration 200-400 ms. (B) ½AM P in response to these bursts was recorded at steady-state, after 4 100 s of simulation. For 4 Hz stimulation, the shorter the burst, the greater the contractile response.
followed by vasodilation. This data indicates that periodically occurring bursts of mSNA dynamically control vascular conductance on a (heart) beat-by-beat timescale, implying that the vasculature is relatively responsive to changes in mSNA. Our model adheres to these data, as it was only able to produce oscillatory contractions in response to bursting input when f r 0:3 Hz (Fig. 7B 1 ). However, our simulation data suggests that the ability of the SNS to transmit oscillatory signals to SMCs depends not only on their incoming frequency, but also their magnitude and duration.

Sympathetic bursting and disease
The results of the present study provide insight into the pathology of hypertension. In the rat, the entrainment of the sympathetic outflow to respiration manifests as bursting activity at respiratory frequencies of 0.5 Hz in situ (Simms et al., 2009;Stalbovskiy et al., 2014). Simms et al. (2009) have demonstrated that the phrenictriggered average of SNA is doubled in amplitude and decreased in duration in the spontaneously hypertensive rat. The sympatheticrespiratory coupling is therefore amplified in this strain. Furthermore, in a recent study , we reported that sympathetic preganglionic neurones in the SH rat display a 40% increase in their average firing frequency (from approximately 2.5 Hz to 3.5 Hz) and an increase of 2-3 spikes per respiratory burst (compared to normotensive controls). These dysfunctions increase respiratory modulation of sympathetic outflow, findings that we have also attributed to increased excitatory respiratory drive to pre-sympathetic neurones in the medulla (Moraes et al., 2014).
In our model, the spike trains from the SH rat were subthreshold for inducing a contractile response, but the respiratory bursts increase NA release and SMC tone. The SH spike trains evoked a 10-fold increase in the contractile response, compared to that evoked by the normotensive train. The amplified respiratory bursts seen in this disease strain would therefore greatly increase SMC tone and contribute to the pathology of high blood pressure.

Model limitations
Our model makes many assumptions about the sympathetic transmission to SMCs, implicit in the theoretical framework. One assumption is that we have played in preganglionic spike trains into a postganglionic model; implicit in this assumption is that there is no pre-to-postganglionic gain. Gain at the level of the sympathetic postganglionic neurones is known to occur due to high-frequency summation of secondary synaptic inputs (Karila and Horn, 2000;Wheeler et al., 2004;Rimmer and Horn, 2010). Given the dependence of secondary gain on high-frequency coincidence and the entrainment of preganglionic activity to respiration, it is reasonable to assume that any gain would occur during the respiratory-modulated burst. Given that the contraction of our SMC model is highly sensitive to changes in burst amplitude (Fig. 7C 1 ), any such gain will markedly influence the contractile response of the cell. The function of postganglionic gain, therefore, could be to add spikes at the critical moment-during bursts -in order to generate a meaningful response from the SMCs comprising the artery wall.
Our model only considers the contractile response due to NA release. Co-transmission is known to occur at sympathetic terminals, with the co-release of neuropeptide Y (NPY) and adenosine 5 0 -triphosphate (ATP) (Pablo Huidobro-Toro and Veronica Donoso, 2004;Burnstock, 2009). Early during a train of postganglionic action potentials, smooth muscle cell contraction is activated by ATP (Wier et al., 2009). ATP is then eliminated from the receptor area in 50-100 ms (Bao et al., 1993). ATP therefore has a short-lived influence on contraction. Later during a train, contraction is maintained due to the binding of NA to α 1 Rs (Wier et al., 2009). Therefore, our model of NA release is likely to provide a good approximation of steady-state contraction of a single SMC to prolonged sympathetic stimulation.
Is the 3.8% increase in SMC tone in response to SH spike trains vs normotensive spike trains biologically significant? This question relates to another aspect of our model-we have modelled the contractile response of a single SMC innervated by a single postganglionic axon. This simplification may explain why the generated contractions from experimental spike trains appear low (o4% of the maximum force available). Considering that SMCs are circumferentially arranged in arteries (Peters et al., 1983;Luff, 1996), a decrease in SMC length translates to a change in arterial radius. As vascular resistance is inversely proportional to the 4th power of this radius (by the Hagen-Poiseuille equation), the small change in SMC length, as seen in response to SH spike trains, can substantially increase vascular resistance. Our future work will model a network of SMCs comprising the arterial wall, innervated by a plexus of sympathetic fibres.

Conclusions
The model has provided a theoretical framework of SMC contraction following stimulation of a sympathetic postganglionic neurone. It provides insight into the dependence of vasoconstriction in arteries on the firing pattern of the neurone. Employment of this model into a macro model of an artery could be used to investigate the influence of altered sympathetic patterning-such as those seen in disease stateson blood pressure. The model predicts that contractility is highly dependent on bursting properties, and that this is due to prejunctional regulation of NA release. Experimental testing of these predictions is important, as they have implications for the altered sympathetic bursting in disease states, such as hypertension (Simms et al., 2009) and heart failure (Sverrisdottir et al., 2000).