Enrico Zio

Roberta Piccinelli

Randomized flow model and centrality measure for electrical power transmission network analysis

Commonly used centrality measures identify the most important elements in networks of components, based on the assumption that flow occurs in the network only along the shortest paths. This is not so in real networks, where different operational rules drive the flow. For this reason, a different model of flow in a network is here considered: rather than along shortest paths only, it is assumed that contributions come essentially from all paths between nodes, as simulated by random walks. Centrality measures can then be coherently defined. An example of application to an electrical power transmission system is presented.

Introduction

Modern society is witnessing a continuous growth in the complexity of the infrastructure networks which it relies upon. Reliable electric power supply, for example, is crucial for many of the services that are taken for granted today; disturbances in the power supply have the potential of severely disrupting these indispensable services. This raises significant concern about reliability and resilience to disturbances and failures of various types of infrastructure systems, and a corresponding demand for methods capable of analyzing the vulnerabilities of these systems [START_REF] Albert | Structural vulnerability of the North American power grid[END_REF]. The developments contained in this paper are motivated by the interest in the analysis of electric power networks. In this context, the network analysis paradigm set up to study the dynamics of the relations in social networks has been previously utilized to analyze the vulnerability of electric power infrastructure systems [START_REF] Holmgren | Using Graph Models to Analyze the Vulnerability of Electric Power Networks[END_REF][START_REF] Zio | Vulnerability analysis of a power transmission network[END_REF]. The focus of these types of studies is typically on analyzing the structural properties of the system from a topological point of view, i.e., considering only the connectivity properties of the network and not the actual physical flow through it [START_REF] Cadini | Using centrality measures to rank the importance of the components of a complex network infrastructure[END_REF][START_REF] Hines | A Centrality Measure for Electrical Networks[END_REF]. Three drawbacks associated with the related measures of network performance are that they are based on:

binary links among network nodes (or components), thus neglecting the strength of the connections (or links or arcs or edges); this has been pointed at as a limitation both in social networks, where the strength and depth of interpersonal relationships is of relevance [START_REF] Festinger | The analysis of sociograms using matrix algebra[END_REF][START_REF] Yan | A fuzzy set analysis of sociometric structure[END_REF][START_REF] Freeman | Centrality in valued graph: a measure of betweenness based on network flow[END_REF] and in engineered network infrastructures, where the capacities of the arcs connecting the components limit the flow among them [START_REF] Hines | A Centrality Measure for Electrical Networks[END_REF]; 2 a simplified modeling scheme which assumes that flow (communication, in the social case) between a pair of components (persons, in the social case) in the network takes place only along the shortest path linking them [START_REF] Freeman | Centrality in valued graph: a measure of betweenness based on network flow[END_REF][START_REF] Stephenson | Rethinking centrality: Methods and examples[END_REF]; this has been considered a limitation in many cases, because the flow from one node of a network to another is typically a global phenomenon which does not depend only on the links on the direct and shortest paths, since it is quite possible that information will take a more circuitous route; this is true both in social networks, where information may travel by random communication or be intentionally channeled through intermediaries, and in network infrastructures, where flow is channeled through selected routes, following the specific operative rules and constraints which apply to the system; a simplified modeling scheme which neglects the possibility of failures in the interconnections between pairs of linked components; this is particularly relevant for the engineered infrastructure networks made of fallible hardware and software, operated by (unfortunately) not error-free human operators.

In synthesis, when looking at the safety, reliability and vulnerability characteristics of an infrastructure such as the electric power transmission network, one should take into account the capacities of the transmission elements and their probability of failure, and examine the different transmission routes available to the flow. This would entail undertaking a complex and detailed mechanistic modeling effort of the entire network system, which is in practice often unfeasible, both with respect to its development and its computation. For this reason, a framework of analysis has been proposed to integrate models at different levels of detail, in a problem-driven approach to solution; complementation of network analysis, for performing an initial screening of the vulnerabilities of a critical infrastructure with object-oriented modeling, to further deepen the vulnerability assessment of the screened scenarios has been investigated as a feasible way to proceed in such direction [START_REF] Eusgeld | The role of network theory and object-oriented modeling within a framework for the vulnerability analysis of critical infrastructures[END_REF]. To improve the physical description of the network characteristics within a network analysis for preliminary screening, a model based on random walks is here introduced as an extension of the model in [START_REF] Freeman | Centrality in valued graph: a measure of betweenness based on network flow[END_REF] giving proper consideration to the following facts:

each link connecting two nodes is characterized by a transmission capacity which cannot be exceeded; the capacities of the network lines are assumed to stochastically vary, to account for the inherent uncertainties; not only the links on the direct and shortest paths are considered in the analysis of the transmission of flow; this is achieved by a randomization of the direction of the flow in output from a node; the randomization is driven by the capacities of the outgoing links, with the highest capacity links most probably channeling the flow; the network interconnecting links are assumed fallible, with given probabilities; source generation and load demands are assumed to vary stochastically, to account for the fluctuations inherent in the network behavior and operation.

From the analysis of the network characteristics and behavior, it is also important to gain an understanding of the role that the elements of the infrastructure network play in determining the flow through it, as this can be of great practical aid to network designers and operators in providing indications for network protection. From a topological viewpoint, various measures of the importance of a network node, can be introduced. These so-called centrality measures, take into account the different ways in which a node interacts/communicates with the rest of the network. Classical topological centrality measures are the degree centrality [START_REF] Nieminen | On Centrality in a Graph[END_REF]12], the closeness centrality [12; 13; 14], the betweenness centrality [12] and the information centrality [START_REF] Latora | A Measure of Centrality Based on the Network Efficiency[END_REF]. The major drawback of these measures is that to assess the node importance they rely only on topological information based on the three previously mentioned model simplifications. Then, based on the model proposed in this paper an extension of the betweenness centrality measure of [START_REF] Ford | Maximal flow through a network[END_REF] is computed, to more realistically capture the importance of the role played by the different components in determining the flow through the network.

An application of the proposed approach is illustrated with reference to a power transmission network system of literature [17]. The paper is organized as follows. In Section 2, a description of the random walk flow propagation model is provided. In Section 3, the topological concept of betweenness centrality measure is recalled and then extended to its randomized flow definition. The results obtained on a case study of literature are discussed in Section 4. Conclusions are drawn in Section 5.

Randomized flow model of a power transmission network infrastructure

The topological interconnection of a power transmission system can be modeled as a network consisting of N nodes (also called vertexes) and K edges (also called arcs or lines): the buses of the electric grid are represented as nodes interconnected by undirected edges representing the transmission lines; N S nodes are power sources (generators), N T nodes are targets (loads) and the rest are transmission nodes. The N×N adjacency matrix {a ij } defines the topological structure of the network, i.e., the pattern of connectivity among its nodes, with the matrix entry a ij being equal to 1 if there is an edge linking i and j and 0 otherwise; the entries on the diagonal elements, a ii , are undefined and for convenience they are set equal to 0. The matrix {q ij } defines the probabilities of failure of the links. The capacities of the links are assumed to vary stochastically, to account for the uncertainties inherent in their behavior and operation; then, to each capacity value w ij is associated a probability distribution 1(w ij ) of the possible values. The underlying strategy to model the flow in the network is to choose a source node, follow at random one of the departing links to one of its neighbors, take this as the source and iterate this process until the required target is reached. The random choice of the arc to follow is based on the actual capacity of each arc outgoing from the node: higher capacity arcs have larger probability to be selected as flow carriers. Accordingly, the algorithm to evaluate the service reliability performance characteristics of the network, and its related vulnerabilities, consists of three nested cycles of randomization; the steps are as follows: 1. Sample the fault configuration of the network on the basis of the failure probabilities of each element (node or arc) of the system. 2. Sample the production from the sources, the demand at the targets and the capacity of the arcs. 3. Build the discrete cumulative distribution function of the capacities of the arcs leaving the source node and sample the flow direction from it. 4. Develop the flow propagation cycle, for each source:

4.1 the random walk of flow follows the arc sampled on the basis of the actual capacities of the arcs departing from the successive nodes traversed by the flow; 4.2 if the flow goes into an isolated node with no departing connections, the cycle ends; 4.3 the flow between a pair of nodes is accounted once (repeated flows between the same pair of nodes are neglected); 4.4 once the flow arrives at a target node, the capacities of the incoming arcs are checked: if their sum is larger than the maximum capacity of the node, an overload is recorded; 4.5 if the flow does not reach the target, a new source of random walk is sampled. If no flow arrives at any of the targets, then a blackout is recorded.

Randomized betweenness centrality measure

Determining the critical elements of large-scale network infrastructures is an important issue for the reliability and the protection of the network. From a topological point of view, a number of centrality indices have been introduced as measures of the importance of the nodes in a network [START_REF] Zio | From Complexity Science to Reliability Efficiency: A new Way of Looking at Complex Network Systems and Critical Infrastructures[END_REF]. These indices take into account the different ways in which a node interacts and communicates with the rest of the network and have proved of value in the analysis and understanding of the role played by the elements in the network. A classical topological centrality measure is the betweenness centrality [12]. This measure is based on the idea that a node is central if it lies between many other nodes, in the sense that it is traversed by many of the shortest paths connecting pairs of nodes. The topological betweenness centrality C i B of a given node i in a network G(N, K), where N is the number of nodes and K is the number of links connecting them, is quantitatively defined as:

1 2 3 4 5 67859678A9 B C DE 629 C DE FF2 (3.1)
where n jk is the number of topological shortest paths between nodes j and k, and n jk (i) is the number of topological shortest paths between nodes j and k which contain node i. C i B assumes values between 0 and 1 and reaches its maximum when node i falls on all geodesics (paths of minimal length between pairs of nodes).

From the definition, betweenness centrality can be regarded as a measure of the influence a node has on the spread of the flow through the network, of the extent to which a node has control over the flow between other nodes. In a network in which flow is entirely or at least mostly distributed along geodesic paths, the betweenness of a vertex measures how much flow will pass through that particular vertex. In most networks, however, flow does not occur only along geodesic paths; in some cases, flow may not follow the ideal route to get from one place to another, and "wander around" in a random-like fashion or as directed by the system operative rules and constraints. In most cases, a realistic betweenness measure should include nongeodesic paths in addition to geodesic ones [START_REF] Borgatti | Centrality and Network Flow[END_REF].

To account for this issue, a betweenness centrality measure based on the concept of network flow has been suggested [START_REF] Freeman | Centrality in valued graph: a measure of betweenness based on network flow[END_REF]. The edges of a network are considered as channels of communication linking pairs of nodes; the value of the connection of two nodes i and j determines the capacity w ij of the channel linking them, or the amount of information that can pass between them. Information is assumed to flow along these channels. If f ij is the amount of information passing on a channel linking node i directly to node j then f ij w ij , i.e., the amount of information flowing along a channel that links directly connected vertices cannot exceed the capacity of that channel. What is relevant here is not just the direct flow between connected nodes, but the overall flow between pairs of nodes along all the paths that connect them: if a node i is chosen as an information source, or transmitter, and another node j as an information sink, or receiver, information from i may reach j along an edge linking i directly to j or along any and all indirect paths that begin at i, pass through one or more intermediate nodes and end at j. Thus, the flow between two nodes is a global phenomenon: it depends, not just on the capacity of the channel linking two nodes directly, but on the capacities of all the channels on all the paths -both direct and indirectthat connect the two. Ford and Fulkerson [START_REF] Ford | Maximal flow through a network[END_REF] introduced a model to determine the maximum flow from any source i to any sink j. Let m jk be the maximum flow from a node i to a node k and let m jk (i) be the maximum flow from j to k that passes through node i; then, the flow betweenness centrality measure may be quantitatively defined as [START_REF] Freeman | Centrality in valued graph: a measure of betweenness based on network flow[END_REF]:

62 2 1 2 3 4 B B E 629 ED D B B DE ED D (3.2)
Flow betweenness is based on the idea of maximum flow; the flow betweenness of a node i is defined as the amount of flow through i when the maximum flow is transmitted from source (s) to target (t), averaged over all s and t. Each edge in a network can be thought of as a transmission line carrying a flow of current. In general, more than a single unit of current can be carried between s and t by making simultaneous use of several different paths through the network. In practical terms, the flow betweenness measures the betweenness of nodes in a network in which a maximal amount of information is continuously pumped between all sources and targets. Necessarily that flow still needs to "know" the ideal route (or one of the ideal routes) from each source to each target in order to realize the maximum flow. This still seems an unrealistic definition, in that it is often the case that flow does not take any sort of ideal path from source to target. To model this, a new betweenness centrality measure has then been introduced, the random walk betweenness [START_REF] Newman | A measure of betweenness centrality based on random walks[END_REF]. Roughly speaking, the random walk betweenness of a node i is equal to the number of times that a random walk starting at s and ending at t passes through i along the way, averaged over all s and t. This measure is appropriate to a network in which information wanders about essentially at random until it finds its target, and it includes contributions from many paths that are not optimal in any sense. Let I i st be the current flowing from s to t, through node i. Quantitatively the random betweenness centrality measure is defined as:

1 2 3 4 B ! " 767859 (3.3)
This measure seems an intuitively reasonable one to describe the fact that current will flow along all paths from source to target, and nodes that lie on no path from source to target get a betweenness of zero. In this paper, an attempt to evaluate random betweenness centrality measures is made considering the physical characteristics of the transmission network in terms of length, capacity and failure probability of each transmission line, and types of nodes.

Application

The artificial transmission network system IEEE 14 BUS [17] is taken as reference case study. The network represents a portion of the American Electric Power System and consist of 14 bus locations connected by 20 lines and transformers, as shown in Figure 1. The transmission lines operate at two different voltage levels, 132 kV and 230 kV. The system working at 230 kV is represented in the upper half of Figure 1, with 230/132 kV tie stations at Buses 4, 5 and 7. The system is also provided with voltage corrective devices in correspondence of Buses 3, 6 and 8 (synchronous condensers). Buses 1 and 2 are the generating units. The failure probability of edge ij is defined as: 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 92 2

# 2F 4 $ % & 8' D ( (4.1)
where 2 ij is the constant failure rate per unit time of the edge ij linking nodes i and j (column 4 in Table II) and T is the reference time for the analysis, here chosen equal to 1 year. Because the failure rate data are usually given as function of the length of each transmission line (column 3 in Table II), in order to compute the failure probability (4.1) transmission line lengths must be known. In this work, these have been inferred from the available data as follows. The total impedance Z ij of a transmission line ij is dependent on the length of the line l ij :

) 2F 4 *+ 2F , -. 2F /0 2F (4.2)
where r ij is the resistance per unit length of arc ij and x ij is its reactance per unit length. While the resistance of the line depends both on the length and on the thickness of the wire, the reactance of the line depends only on the length [START_REF] Grisby | Power Systems[END_REF]. Based on relation (4.2), the lengths of the transmission lines in the IEEE 14 BUS system have been obtained from literature data by taking a system power base of 100 MVA and a conversion factor of 0.48 3/km [START_REF][END_REF]; lines containing transformers are considered to be zero-length lines (Table III).

From BUS i

To The shortest lines are concentrated in the lower half of the network, which contains the generating units, while the longest lines belong to the upper half of the network, which contains the loads. The largest failure probabilities are concentrated on the edges directly connected with sources and with loads (edges 1 -5, 2 -3, 2 -4 and 12 -13, 9 -14, 13 -14).

A2 2

Source generation is sampled from a normal distribution with a mean value of 30 and a variance of 100, in arbitrary units (a.u.). The values of the capacities of the network links are assumed all distributed according to a normal distribution of mean value 100 a.u. and a standard deviation of 10 a.u. The direction of flow is sampled on the actual capacities of the arcs. Once the flow arrives at a target node, the capacities of the incoming arcs are checked: if their sum is larger than the maximum capacity value of the node, an overload is recorded. The targets are absorbing nodes: the flow stops and the received flow is recorded for evaluating the network lost load and the network service efficiency; if no flow reaches any target, a service blackout is recorded. The network performance characteristics computed on the basis of the above data are reported in Table IV where:

-blackouts and overloads are evaluated considering the average value of the flow that does not reach the targets or that exceeds the capacities of the transmission lines, respectively; -the network demanded load is the average sum of the power generated from all the sources s i , i=1, 2,…, N S :

123 4 B 4 2 7 5 265 (4.3)
-the network received load is the average sum of the flow reaching the targets t i , i=1, 2, …, N T :

13 4 B 7 2 7 8 265 (4.4)
-the network lost load is obtained as the difference between demanded and received loads:

133 4 123 % 13 (4.5)
-the network service efficiency is obtained as the ratio between received and demanded loads:

19: 4 7;< 7=< (4.6) 1B2 2
In the artificial case considered, the network service efficiency is low and the blackout probability is high: very little of the generated power is received from the loads. This is caused by the high probability of failure of the transmission lines, due to their relatively high failure rates and large lengths. To see the effect of the lines lengths, a second computation has been made using average transmission line lengths inferred from literature [START_REF][END_REF]. Two line lengths of 48 and 50 km have been considered (Table V With respect to the previous case, targets now receive a larger load, the service efficiency is higher and the blackout probability is lower. This is due to the values of the transmission lines failure probabilities, which are smaller than in the previous case because of the smaller line lengths. Finally, the centrality betweenness based on the proposed random walk model has been computed for the two examples. The results are shown in Figure 3. Betweenness values evaluated for case 1 (×, Figure 3) using higher lengths computed from eq. (4.2), are equal or higher than the values obtained for case 2 using smaller average lengths (5, Figure 3), due to the associated larger failure probabilities. Equal values are obtained for nodes 1 and 2, i.e., the source nodes, for node 3 and for node 8 which is an isolated node (Figure 1). The lower half of the network (nodes 1, 2, 3, 4 and 5), which contains the generating units (nodes 1 and 2), has higher values of betweenness than the upper half, which contains the load nodes 13 and 14. Nodes 10, 11 and 12 act as a tie for the flow. Node 2, which represents a source node and a transfer node as well, is the most important from the betweenness point of view, in both studies.

Conclusions

In previous works, network analysis has been shown suitable for a preliminary analysis of complex infrastructures aimed at identifying structural criticalities, e.g. the most connected nodes, shortest path lengths of connection, most vulnerable nodes, etc. Limitations of the analysis relate to the neglecting of the actual capacities of the links, their probabilities of failures and the fact that flow among network nodes is typically a global phenomenon, not restricted to only direct, shortest paths as typically assumed. To overcome some of these limitations, in this paper a model of random flow propagation has been introduced and the topological concept of betweenness centrality has been accordingly extended to account for the random flow propagation across the network. The randomization of the flow out of a node is driven by the capacity values of its outgoing links and allows non-geodesic paths to be travelled by the flow. Variability in the behavior and operation of the links, source and target nodes is also accounted for by varying stochastically the capacities, the productions and demands, respectively. The modeling approach has been applied to the artificial transmission network system of the IEEE 14 BUS and indications derived from the betweenness centrality measure values have been analyzed for different transmission lines lengths. Each equipment of the system has been transposed into a node or edge of the representative network and the length of the arcs has been calculated for two cases: case 1, in which the lengths have been obtained from the impedances of the lines and case 2, in which average line lengths have been considered from literature data. The network performance characteristics and the random walk betweenness centrality measures have highlighted the weaknesses of the network structure, for the failure data used.

Figure 3 .

 3 Figure 3. Betweenness centrality according to the random flow model proposed.

Table I .

 I the adjacency matrix that defines the topological structure of the network. The IEEE 14 BUS adjacency matrix Table II reports the failure rate values of the components of the transmission network, as inferred from literature data[START_REF] Billington | Reliability Assessment of electric power system using Monte Carlo Methods[END_REF].

			{a ij } =
	1
	1
	1
				Fig Figure 1. Transm smission network rk IEEE 14 BUS S [17].
	To carry ou out the analys lysis, each net network comp ponent is tran transposed into nto a node or or edge of the the representat tative
	topological al network, as as it is shown n in Figure 2. 2. Three diffe fferent physica ical types of no nodes are con onsidered: sou source
	nodes (wher here the electric tricity is fed in into the netwo work), load nod nodes (where c e customers are are connected) ed) and transfer sfer or
	transmission ion nodes (with ithout custome mers or source rce).
	2
	From BUS To BUS Failure rate/yr/100km Failure rate/yr	Equipment
	1		2	1.0858	132 kV transmission line
	1		5	1.0858	132 kV transmission line
	2		3	1.0858	132 kV transmission line
	2		4	1.0858	132 kV transmission line
	2		5	1.0858	132 kV transmission line
	3		4	1.0858	132 kV transmission line
	4		5	1.0858	132 kV transmission line
	4		7	0.01045	132/230 kV transformer
	4		9	0.01045	132/230 kV transformer
	5		6	0.01045	132/230 kV transformer
	6		11	0.5429	230 kV transmission line
	6		12	0.5429	230 kV transmission line
	6		13	0.5429	230 kV transmission line
	7		8	0.01045	132/230 kV transformer
	7		9	0.01045	132/230 kV transformer
	9		10	0.5429	230 kV transmission line
	9		14	0.5429	230 kV transmission line
	10		11	0.5429	230 kV transmission line
	12		13	0.5429	230 kV transmission line
	13	1	14	0.5429	1 230 kV transmission line
	Figure 2. The he IEEE 14 BUS S transmission n n network's graph aph representation ion [21]. The whi hite circles label beled with G repr epresent the gener
	nodes, the colo olored circle node des represent the he transmission n n nodes and the w white cylinders l s labeled with L	represent the lo load nodes.
				72

nerator Table

I provides

Table II .

 II Failure data of the arcs of the IEEE 14 BUS transmission network

Table III .

 III Length and failure probabilities of the transmission lines

		BUS	Length (km)	Failure Probability
		j	l ij	q ij
	1	2	22	0.2125
	1	5	81	0.5768
	2	3	72	0.5338
	2	4	64	0.4932
	2	5	63	0.4884
	3	4	62	0.4828
	4	5	15	0.1498
	4	7	-	0.0104
	4	9	-	0.0104
	5	6	-	0.0104
	6	11	220	0.6970
	6	12	283	0.7847
	6	13	144	0.5425
	7	8	-	0.0104
	7	9	-	0.0104
	9	10	93	0.3978
	9	14	299	0.8027
	10	11	212	0.6843
	12	13	221	0.6999
	13	14	385	0.8762

Table IV .

 IV Network performance indicators

		:
	Blackout (%)	0.82
	Overload (%)	0.00
	Network service efficiency	0.18
	Network demanded load (a.u.) 59.93
	Network received load (a.u.)	10.92
	Network lost load (a.u.)	49.01

Table V .

 V ). Length and failure probabilities of the arcs In this case, the lower half of the network has the largest failure probability. The network performance characteristics are reported in TableVI.

	From BUS	To BUS	Length (km)	Failure probability
	i	j	l ij		q ij
	1	2	48		0.4079
	1	5	48		0.4079
	2	3	48		0.4079
	2	4	48		0.4079
	2	5	48		0.4079
	3	4	48		0.4079
	4	5	48		0.4079
	4	7	-		0.0104
	4	9	-		0.0104
	5	6	-		0.0104
	6	11	50		0.2372
	6	12	50		0.2372
	6	13	50		0.2372
	7	8	-		0.0104
	7	9	-		0.0104
	9	10	50		0.2372
	9	14	50		0.2372
	10	11	50		0.2372
	12	13	50		0.2372
	13	14	50		0.2372
	Blackout (%)			0.44
	Overload (%)			3.33410 -4
	Network service efficiency		0.60
	Network demanded load (a.u.)	59.93
	Network received load (a.u)		36.24
	Network lost load (a.u)		23.70

Table VI .

 VI Network performance indicators

	112
	2

Acknowledgements

This work has been partially funded by the Foundation pour une Culture de Securité Industrielle of Toulouse, France, under the research contract AO2006-01.

Acknowledgments:

The authors wish to thank Prof. Maurizio Delfanti and Dr.Mauro Pozzi of the Department of Energy of the Politecnico di Milano, for their contribution to the work.


Enrico Zio
email: enrico.zio@polimi.it

Francesco Di Maio

A FUZZY SIMILARITY-BASED METHOD FOR FAILURE DETECTION AND RECOVERY TIME ESTIMATION

Keywords: Failure Detection, Recovery Time Estimation, Emergency Accident Management, Lead-Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS), Fuzzy Similarity Analysis, Pattern Recognition

In this work, an extension of a data-driven approach for estimation of the available Recovery Time (RT) is presented. The improvement is in avoiding the need of resorting to a fault detection module for the identification of the anomalous system behavior: the algorithm proposed jointly detects the onset of the accidental transient and performs the estimation of the available RT. This is achieved by fuzzy similarity analysis of the currently developing scenario and reference multidimensional trajectory patterns of failure scenarios; the RT remaining before the developing trajectory pattern hits a failure threshold is predicted by combining the times of failure of the reference patterns, weighed by their similarity with the developing pattern.

For illustration purposes, failure scenarios of the Lead Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) are considered.

Introduction

Nuclear Power Plant (NPP) accident management involves the anticipation of paths of potentially dangerous behaviors, the prediction of the related effects and the setting of actions to avoid undesired impacts on the safety of the NPP [IAEA, 2003].

In case of an accident, or an initiating event that may develop into an accident, the plant personnel must perform various tasks before taking counteracting actions:

• Identification of the plant state: this diagnostic task aims at identifying the cause of the problem and the states of a number of parameters critical for the plant operation and safety;

• Prediction of the future development of the accident: this involves prediction of the future evolution of the states of the critical parameters and of the plant residual Recovery Time (RT), which is the time remaining for taking corrective actions to avoid system failure;

• Planning of accident mitigation strategies, to be activated if the actions to safely control the accident evolution were not successful.

Computerized procedures can be used to aid the operators in the tasks of obtaining reliable information, following procedures, identifying plant states, predicting the future accident progression and planning accident control and mitigation actions. Indeed, it is a shared belief that the complicated phenomena taking place in a NPP during an accident situation need to be supported by computerized procedures rather than left only to human expert evaluation [Øwre, 2001]. On the other hand, the problem of what kind of support to provide to nuclear power plant operators, in particular during transients leading up to accidents, is far from trivial [NEA, 1992;[START_REF] Ec | EU research in Reactor safety: Achievements of the 4th and prospects for the 5Th Euratom framework programme[END_REF]USNRC, 1999;IAEA, 2003].

Fortunately, NPP personnel have the capability to effectively manage a broad range of accidents; their successful handling of complex accident behaviors requires that they detect the occurrence of the accident, determine the extent of challenge to plant safety, monitor the performance of active, passive, automatic and digital systems, select strategies to prevent or mitigate the safety challenge, implement the action strategies, and monitor their effectiveness. The capability to effectively carry out these tasks during an accident is influenced by the availability of timely and accurate plant status information and the awareness of the RT available for action. Poor decisions may be taken under a perceived short time available for sorting out the information relevant to the plant status [START_REF] Glasstone | Nuclear reactor engineering[END_REF]; on the contrary, timely and correct decisions can prevent an event from developing into a severe accident or mitigate its undesired consequences.

The existing computer-based tools which can aid accident management can be categorized according to their complexity and specific application purpose [IAEA, 2003]:

• Compact simulators;
• Plant analyzers;

• Full-scope training simulators;

• Multifunctional simulators;

• Severe accident simulators;

• Accident management support tools.

In particular, the accident management support tools combine tracking and predictive simulators of the developing accidental transients. The tracking simulators monitor the plant status and provide calculated values also of those parameters that are not directly observable by the monitoring systems.

One of the main quantities of interest to be delivered by predictive simulators for accident management use is the available RT. Approaches to RT prediction can be categorized broadly into model-based and data-driven [START_REF] Chiang | Fault detection and diagnosis in industrial systems[END_REF]. Model-based approaches attempt to incorporate physical models of the system into the estimation of the RT. However, uncertainty due to the assumptions and simplifications of the adopted models may pose significant limitations. Moreover, the requirement of high computational speed for on line response necessarily leads to limited details in the phenomena modeled, with consequent limited accuracy in the representation of the actual plant behavior [Berglund et al., 1995;Serrano et al., 1999]. On the contrary, data-driven techniques utilize monitored operational data related to system health. They can be beneficial when understanding of first principles of system operation is not straightforward or when the system is so complex that developing an accurate model is prohibitively expensive. Data-driven approaches can often be deployed quickly and cheaply, and still provide wide coverage of system behavior.

Data-driven approaches can be divided into two categories: statistical techniques (regression methods, ARMA models, etc.) and Artificial Intelligence (AI) techniques (neural networks, fuzzy logic systems, etc.).

With respect to AI techniques, Neural Networks (NNs) and Fuzzy Logic (FL) techniques have gained considerable attention in the past few years, due to their ability to deal with the uncertainties and non-linearities of the real processes, especially in abnormal conditions [Øwre, 2001]. Successful prediction models have been constructed based on Neural Networks [Barlett et al., 1992;[START_REF] Campolucci | On-Line Learning Algorithms of Locally Recurrent Neural Networks[END_REF]Peel et al., 2008;[START_REF] Zio | Locally recurrent neural networks for nuclear dynamics modeling[END_REF][START_REF] Santosh | Diagnostic System for Identification of Accident Scenarios in Nuclear Power Plants using Artificial Neural Networks[END_REF] and Neuro-Fuzzy (NF) systems [START_REF] Wang | Prognosis of Machine Health Condition using Neuro-Fuzzy Systems[END_REF]. In spite of the recognized power of neural network modeling techniques, skepticism against their use in safety-critical applications relates to their black-box character which limits intuition with respect to the understanding of their performance [START_REF] Wang | A similarity based prognostic approach for Remaining Useful life estimation of Engineered Systems[END_REF]].

An opportunity for increased transparency and openness of data-driven models is offered by fuzzy logic methods, which are increasingly proposed in modern control and diagnostic technologies. Based on the principles of Zadeh's fuzzy set theory, fuzzy logic provides a formal mathematical framework for dealing with the vagueness of everyday reasoning [Zadeh, 1965]. As opposed to binary reasoning based on ordinary set theory, measurement uncertainty and estimation imprecision can be properly accommodated within the fuzzy logic framework [START_REF] Yuan | Data driven identification of key variables[END_REF][START_REF] Zio | Identification of nuclear transients via optimized fuzzy clustering[END_REF].

The work presented in this paper extends the multidimensional approach previously presented by the authors for the prediction of the available RT during an accident [Zio et al., 2009b]. In particular, the need is avoided of a fault detection module that from the beginning matches the evolving trajectory with the trajectories stored in the database; in view of a simplification of the procedure for its practical application and an increase of the reliability of the overall accident management support system, the RT estimation module is devised so as to automatically identify the anomalous behavior of the plant and proceed to the RT evaluation.

The computational framework considers a set of multidimensional trajectory patterns arising from different system failures (hereafter called reference trajectory patterns) and uses a fuzzybased, data-driven similarity analysis for predicting the RT of a newly developing failure trajectory (hereafter called test trajectory pattern). The pattern matching process is based on a fuzzy evaluation of the distance between the signals of the multidimensional test pattern and the patterns of reference [Angstenberger, 2001]; the fuzzy distances from all reference patterns are then combined to transform the multidimensional data into a one-dimensional similarity indicator, which is used in the estimation of the available RT.

An application of the method is illustrated with regards to the dynamic failure scenarios of the Lead Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) with digital Instrumentation and Control (I&C) [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF]. The focus of the application is on the RT estimation and, therefore, the modeling is tailored to the purpose of showing the feasibility of effectively identifying the onset of a failure scenario and estimating the RT during an accident: thus, the analysis does not cover the software failure behavior, the benefits of fault tolerant features, the interactions of the software with the hardware and human components; it is understood that the actual implementation of the method in a qualified tool of accident management will require full dynamic accident calculations, inclusive of comprehensive models of hardware, software and human failure modes and their interactions.

The paper is structured as follows. Section 2 provides a detailed description of the computational algorithm for the multidimensional fuzzy data similarity evaluation and the associated RT prediction. Section 3 presents the mechanistic model of the LBE-XADS, with the description of the monitored signals. In Section 4, the results of the application of the approach to LBE-XADS failure scenarios are presented. Finally, some conclusions are drawn in Section 5.

Methodology of the RT estimation by fuzzy similarity analysis

It is assumed that N trajectories (reference trajectory patterns) are available, representative of the evolution of relevant signals during reference failure scenarios. These trajectories last all the way to system failure, i.e., to the time when anyone of the signals reaches the threshold value beyond which the system loses its functionality.

At any time, a multidimensional signal trajectory (test trajectory pattern) is compared for similarity with the N multidimensional reference trajectory patterns stored in the database and the residual lifes of these are used to estimate the RT available in the developing transient.

Figure 1 shows a schematics of the computational framework in the general case of multidimensional trajectories of Z monitored signals ( ) For completeness of the paper, procedural Steps 1-5 of the algorithm, already presented in [Zio et al., 2009a], are reported in Appendix. Here, only the key modifications with respect to the previously presented procedure (i.e., the addition of Step 6) are highlighted, together with the motivation behind it.

1 2 , ,..., , Z f x x x t .
-Step 1: signal monitoring.

See Appendix.

-Step 2: pointwise difference computation of multidimensional trajectory patterns.

See Appendix.

-Step 3: computation of trajectory pointwise similarity and corresponding distance score.

See Appendix.

-Step 4: weight definition.

See Appendix.

-Step 5: RT i (t) and RT(t) estimation.

See Appendix.

-Step 6: Stationary signal identification and RT(t) estimation as MTTF(t).

Starting from time 1 t = , at each discrete time t , the devised algorithm is qualified to continuously provide the prediction ( ) RT t . However, in case that the developing signal is actually stationary, ( ) RT t would be biased towards too conservative estimations because of the "max" assumption adopted in Eq. ( 10

( ) ( ) ( ) ( ) ( ) ( ) | | 1 1 | | | | i f f i i i i f i i t t i t t f f RT t MTTF t t t RT t i t t i t t > > = = -= > > ∑ ∑ (1)
where i f t is the system failure time along the i-th trajectory, ( ) In the subsequent time steps, as the scenario evolves, the procedure restarts from Step 1 and continues to update the estimate of ( ) RT t : in case that non-stationarity in the signal is detected, the ( ) RT t is computed through steps 1-5 of the procedure; otherwise,

| | i f i t t >

( )

RT t is taken equal to the ( ) MTTF t .

Case study

The case study considered is the same as in [Zio et al., 2009a], concerning the RT in failure trajectories of the Lead-Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS), a sub-critical, fast reactor in which the fission process for providing thermal power ( ) P t is sustained by an external neutron source through spallation reaction by a proton beam ( ) Q t accelerated by a synchrotron on a lead-bismuth eutectic target [START_REF] Bowman | Nuclear energy generation and waste transmutation using an accelerator-driven intense thermal neutron source[END_REF][START_REF] Carminati | An Energy Amplifier for Cleaner and Inexhaustible Nuclear Energy Production Driven by a Particle Beam Accelerator[END_REF][START_REF] Rubbia | Conceptual Design of a Fast Neutron Operated High Power Energy Amplifier[END_REF]Van Tuyle et al., 1993;[START_REF] Venneri | Accelerator-driven Transmutation of Waste (ATW) -A New Method for Reducing the Long-term Radioactivity of Commercial Nuclear Waste[END_REF]. The physical description of the LBE-XADS and of the failures considered are here repeated, for completeness of the information in the paper. The primary cooling system is of pool-type with Lead-Bismuth Eutectic (LBE) liquid metal coolant leaving the top of the core, at full power nominal conditions, at temperature equal to 400 °C and then re-entering the core from the bottom through the down-comer at a temperature of 300 °C.

The average in-core temperature of the LBE , av C LB T is taken as the mean of the entrance and exit temperatures. The secondary cooling system is a flow of an organic diathermic oil at 290-320 °C, at full power conditions. Cooling of the diathermic oil in each loop is obtained through an air flow ( ) a t Γ provided by three air coolers connected in series.

A dedicated, dynamic simulation model has been implemented in SIMULINK for providing a simplified, lumped and zero-dimensional description of the coupled neutronic and thermohydraulic evolution of the system [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF]. The model allows the simulation of the system controlled dynamics as well as of the free dynamics when the control module is deactivated and the air cooler flow is kept constant.

Both feedforward and feedback digital control schemes have been adopted for the operation of the system. The control is set to keep a steady state value of approximately 300 °C of the average temperature of the diathermic oil T =280 °C could result in thermal shocks for the primary fluid and, eventually, for the structural components [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF]. Conservatively, no dependence on the duration of exposition to temperatures beyond the threshold values has been assumed: in other words, the system is considered to fail at such temperatures regardless of the time during which it exceeds the thresholds.

Multiple component failures can occur during the system life. To simulate this, the model has been embedded within a Monte Carlo (MC) sampling procedure for injecting faults at random times and of random magnitudes. Samples of component failures are drawn within a time horizon of 3000 s. The set of faults considered are:

• The PID controller fails stuck, with a random flow rate output value 1 m sampled from a uniform distribution in [0,797] Kg/s.

• The air coolers fail stuck in a random position that provides a corresponding air flow mass 2 m uniformly distributed in [0,1000] Kg/s.

• The feedforward controller fails stuck with a corresponding flow rate value 3 m uniformly distributed in [0,797] Kg/s.

• The communication between air coolers actuators and PID controller fails so that the PID is provided with the same input value of the previous time step.

The first three faults are applicable to both analog and digital systems, whereas the last one is typical of digital systems. Furthermore, the fault magnitude probability distributions are assumed to be uniform, even if the components may more likely fail in a certain mode than in others. This includes also rare multiple events in the set of failure scenarios and further tests the robustness of the RT prediction procedure.

The sequence of multiple failures is generated by sampling the first failure time from the uniform distribution [0,3000] s and the successive failure times from the conditional distributions, uniform from the last sampled time to 3000 s. This assumption is conservative, favoring larger number of failures in the sequence [START_REF] Maio | [END_REF].

The evolution of the failure scenarios may lead to three different end states, within the time horizon of 3000 s:

1. Low-temperature failure mode ( , It is important to stress that the procedure for sampling the fault events here implemented is not intended to reproduce the actual stochastic failure behavior of the system components; rather, the choices and hypotheses for modeling the faults (i.e., the time horizon of the analysis, the number and typology of faults, the distributions of failure times and magnitudes) have been arbitrarily made with the aim of favoring multiple failures. In any case, the components considered subjected to fault and their types are not intended to provide a comprehensive description of the system fault behavior but are only taken as exemplary for generating the dynamic failure scenarios to be used as reference and test patterns.

Results

Application of the procedure for RT(t) prediction

For the present case study, Z=2-dimensional reference trajectories of evolution of T , relative to N = 6400 failure scenarios (differing in faulty components, times of faults occurrence and faults magnitudes), have been used as in [Zio et al., 2009b]. The database of reference trajectories is organized in the structure

[ ] N k Z R × × , where T k n = =30. The generic element ( )
, , r i j z of the reference structure will be compared for similarity with the z-th signal of the test trajectory pattern containing the values of the latest n=100 time steps of the trajectory. For each of the test trajectories the six steps of the procedure (see Appendix) are performed as follows:

Step 1: signal monitoring.

The 2-dimensional trajectory ( )

1 2
, , f x x t is followed starting from time 1 t = s to the end of the time horizon T =3000 s, with time steps of 1 s. At each time step t , its value is appended and stored in the matrix containing the 1 n -=99 values of the two trajectory signals collected at the previous times. The fuzzy similarity algorithm for RT estimation is activated.

Step 2: pointwise difference computation of multidimensional trajectory patterns.

The total pointwise difference ( ) , i j δ is evaluated by resorting to Eq. (3) in the Appendix, with z=1, 2.

Step 3: computation of trajectory pointwise similarity and corresponding distance score.

The pointwise differences ( ) , implying strong sharpness in the FS and thus in the similarity requirement [Zio et al., 2009a]. The distance scores ( ) , d i j are then computed by Eq. ( 4) in the Appendix, 1, 2,..., 6400 i = , 1, 2,..., 30. j =

Step 4: weight definition.

The minimum distances * i d are evaluated (Eq. ( 7) in the Appendix), and the relative normalized weights i w calculated through Eqs. ( 8) and ( 9) in the Appendix, 1, 2,..., 6400. i =

Step 5: RT i (t) and RT(t) estimation.

For each reference trajectory in the library, an estimate ( ) i RT t for the test trajectory is computed (Eq. ( 10) in the Appendix, 1, 2,..., 6400 i =

); then, the values ( ) i RT t are aggregated in the weighted sum (Eq. ( 6) in the Appendix) with the weights i w previously calculated.

Step 6: stationary signal identification and RT(t) estimation as MTTF(t).

In case a stationary signal is detected (i.e., ( ) ( ) ( ) ( )

1 2 RT t RT t RT t RT t l = -= -= -, 1, 2,..., l T =
), the Mean Time to Failure ( ) MTTF t is calculated resorting to Eq. ( 1) and

( )

RT t is set equal to it, for each time step.

The estimates of the ( )

MTTF t for five 2-D test pattern trajectories taken from [Zio et al., 2009a] are plotted in Figures 34567, in thin continuous lines with the bars of one standard deviation of the samples ( )

| i i f f t t t t - > ,

Performance evaluation of the RT(t) estimation procedure

The performance of the available RT estimation procedure has been tested extensively on a batch of P=1280 multidimensional test trajectories, different from the reference trajectories in the fuzzy similarity analysis upon which the RT estimation is based. Figure 8 shows the ( )

RT t predictions (continuous line with dark bullets) compared to the actual remaining lives for 25 failure trajectories (graphically appended in sequence), starting from the time in which the trajectories deviates from the stationary nominal value. The assumption made in step 5 of the procedure for the evaluation of Eq. ( 7) in the Appendix results in a conservative trend of initial anticipation of the available RT associated to trajectories whose failure actually occurs late in life. Also, as expected, the largest available RT estimation errors occur for those trajectories in which the component failure is of low-magnitude, whose effect only slowly drives the system to failure and .

Figure 9 shows the empirical probability density function of the mean relative error evaluated on 3000 s of evolution of the 1280 test trajectories; the distribution is skewed towards small error values, with mean and median equal to 0.08 and 0.03, respectively: this proves that the procedure most frequently makes small relative estimation errors.

The computational time required for the estimation along one complete test trajectory of 3000 s is of 20 s on an Intel ® Celeron M of 900 MHz. [Zio et al., 2009b] while in the other case the results have been obtained with the algorithm as implemented in this work, i.e., without resorting to any "a priori" fault detection module. In both cases, β in Eqs. ( 5) and ( 8) of the Appendix is set equal to 0.02. It is seen that:

-the accuracy in the estimation of the available RT improves over time: as the available RT decreases, the relative errors approximate their mean values approximately equal to 0.05.

-In the latest steps, the accuracy in the estimation of the available RT never exceeds a RE value of 0.1.

-The accuracy of the unattended RT estimation algorithm here implemented is comparable with the one that resorts to the fault detection module for the identification of an anomalous transient.

-Also, it is expected that, by avoiding the fault detection task, the RT estimation procedure is made free of possible errors in the fault detection module, so that the reliability of the whole operator support system for emergency accident management is increased.

Conclusions

This paper extends a fuzzy similarity-analysis method for estimating the available Recovery Time (RT) during the evolution of a failure trajectory of a system.

The extension to the method is aimed at freeing it from the need of resorting to a fault detection module for the identification of the anomalous system behavior. This is expected to increase the reliability of the method by avoiding potential errors in the detection module. The method of RT estimation relies on data from different transient failure scenarios to create a library of reference patterns of evolution; for estimating the available RT of a test pattern, its evolution data are, then, matched to the patterns in the library and their known residual life times are used for the estimation, based on a multidimensional fuzzy pointwise similarity analysis procedure.

The RT estimation procedure involves six main steps: (1) monitoring the signal evolution;

(2) computing the pointwise difference between the test pattern evolution monitored and the available reference patterns;

(3) evaluating their multidimensional fuzzy pointwise similarity and distance score; (4) defining the weights of the individual RT estimates provided by the fuzzy similarity of the reference patterns; (5) aggregating these to evaluate the system residual RT; and

(6) setting to the MTTF value the RT estimation, when the system is identified as working in stationary conditions.

Application to a literature case study concerning the RT estimation for fault scenarios in the Lead Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) has given satisfactory results from the point of view of both accuracy of the estimation and computing time.

The increased computing time required by matching the trajectories at each time step, and not just after the fault detection module has indicated a fault, is expected to be repaid by the increased reliability of the RT estimation system. Reliability and computational time are two objectives which can be compromised case-by-case through a proper selection of the number of monitored signals upon which to base the pointwise similarity evaluation and with the use or not of a fault detection module. At the current time t, the latest n-long segment of values of the test trajectory pattern ( ) ( )

1

-Step 3: computation of trajectory pointwise similarity and corresponding distance score.

In general, there are numerous practical cases in which the classification of objects as 'similar' or 'non similar' is gradual; in these cases, then, the similarity measure should allow for a gradual transition [START_REF] Binaghi | Fuzzy Reasoning Approach to Similarity Evaluation in Image Analysis[END_REF][START_REF] Joentgen | Dynamic Fuzzy Data Analysis based on Similarity between Functions[END_REF]. In the case of interest here, this can be achieved by evaluating the pointwise difference of two trajectories with Note that the smaller the minimum distance, the larger the weight given to the i-th trajectory [Zio et al., 2009a].

-Step 5: RT i (t) and RT(t) estimation.

Referring to the generic i-th trajectory in the library for which the system failure time i f t t > (i.e., the time when the signal value exceeds the threshold beyond which the system loses its functionality), the value

Figure 1

 1 Figure 1 The flowchart of the procedure for RT estimation by fuzzy similarity analysis

  is the cardinality of the set of reference trajectories whose failure time is larger than t and ( ) i RT t is their residual life starting from t .

Figure 2

 2 Figure 2 The flowchart of the procedure for RT estimation by fuzzy similarity analysis applied to a stationary signal

  value represents the optimal working point of the diathermic oil at the steady state, full nominal power of 80 MWth. On the contrary, an oil temperature beyond the upper threshold , th u o T =340 °C would lead to degradation of its physical and chemical properties, whereas a temperature below the lower threshold , th l o

••

  signals are taken for the estimation of the available RT:• Mean in-core LBE temperature, Mean oil temperature of the secondary heat exchanger hot side , Mean air flow rate at the secondary heat exchanger cold side, ( ) a t Γ

  " FS. The bell-shaped function of Eq. (5) in the Appendix is taken

  where i f t is the time at which the diathermic oil temperature profile of the i-th reference trajectory exceeds either thresholds , bold circles. Initially, the estimate of the RT matches the ( ) MTTF t , because the algorithm has not yet acquired enough information about the developing trajectory; then, the algorithm starts to provide the ( ) RT t , based on the results of the fuzzy similarity analysis: the ( ) RT t values matches the ( ) MTTF t in case of a stationary signal, or drifts away from the ( ) MTTF t values and towards the real RT(t) (dashed thick line), if an anomalous signal deviation is detected. In the Figures, the bold vertical line indicates the time of diathermic oil threshold exceedance; notice that none of the estimates exceeds the actual failure time.

Figure 3

 3 Figure 3 RT estimation for a trajectory belonging to the lowtemperature failure mode

FigureFigure 9

 9 Figure 8 Predicted and actual remaining life for 25 test failure scenarios

Figure 10

 10 Figure 10 Relative Error evaluated each 50 s starting from 600 s before failure, for 1280 test trajectories, using Z=2 signals: in one case the fault detection module is activated [Zio et al., 2009b] (circles), in the other case the fault detection module is not implemented (squares)

Figure

  Figure 11 A bidimensional reference trajectory on a time horizon T=3000 s

  n-long segments of the Z-dimensional reference trajectories and the test trajectory pattern of the monitored signals.

…

  trajectory should be more informative in the extrapolation of the occurring trajectory to failure.To assign the weight i w , the minimum distance * i d along the i th row of the matrix of Eq. (4) is first identified:

  time index of the latest-in-life segment of the i-th trajectory among those with minimum distance * i d from the developing test

  ) in the Appendix for the computation of horizon of observation T ) and, therefore, iii) the correction of the bias of extremely low (i.e., pessimistic) RT predictions. With respect to the latter point, when at time t the signal is classified as stationary, the correction consists in the substitution of the estimate

	( ) RT t with the system Mean Time to Failure	( ) MTTF t , obtained from the available
	recovery times	( )
									M j t ; under
	this	assumption,	an	(l+1)-long			stationary	segment	entails	that
	( ) RT t RT t ( ) 1 = -=	( ) 2 RT t -=	( ) RT t l -,	l	=	1, 2,...,	T	, as shown in Figure 2. The
	detection of such string of RT estimates allows for i) the classification of the signal as of no

abnormal deviation ii) the corresponding qualification of the system state as working in nominal conditions (in other words, no fault has been detected and the system reaches safely the time i RT t of the entire population of failure trajectories in the reference library:

Acknowledgements

This work was supported in part by Ordine degli Ingegneri di Milano and by the Science & Technology Fellowship Programme in China financed by the European Commission (Reference: EuropeAid/127024/L/ACT/CN).

APPENDIX

The procedure of fuzzy similarity analysis and RT estimation follows according to six steps:

-Step 1: signal monitoring.

The Z-dimensional trajectory of signals

, ,..., , Z f x x x t is continuously monitored throughout the time horizon of observation T , starting from (discrete) time 1 t = ; at each discrete time t , its values are recorded and appended to the matrix of the values collected at the previous time steps. For reasons which will become clear in the following, the database containing the reference trajectory patterns is organized in a

= and its generic element ( ) , , r i j z is the projection on the z-th signal axis of ( ) , r i j which is the j-th segment, of length n , of values of the i-th reference trajectory, At each time, the algorithm for the estimation of the available RT matches the similarity of the developing signal trajectory to those in the reference library and combines their end times to provide an estimate of the available RT of the current trajectory.

reference to an "approximately zero" fuzzy set (FS) specified by a function which maps the elements ( )

the condition of "approximately zero". The distance score ( ) , d i j between two trajectory segments is then computed as:

Common membership functions can be used for the definition of the FS, e.g. triangular, trapezoidal, and bell-shaped [START_REF] Dubois | Weighted Fuzzy Pattern Matching[END_REF]. In the application illustrated in this work, the following bell-shaped function is used:

The arbitrary parameters α and β can be set by the analyst to shape the desired interpretation of similarity into the fuzzy set: the larger the value of the ratio ( )

narrower the fuzzy set and the stronger the definition of similarity [Zio et al., 2009a].

-Step 4: weight definition.

The ( )

RT t is estimated as a similarity-weighted sum of the

The ideas behind the weighting of the individual Then, the estimate ( )

RT t of the remaining useful life along the developing trajectory is simply computed as in Eq. ( 6), with weights i w evaluated by Eq. ( 9).

It has to be noted that the assumption

allows for a conservative RT estimation, biased towards "pessimistic" predictions of the available RT, because in the case that more than one segment along the i-th reference trajectory is closest to the developing test trajectory, the latest one is taken, i.e., the one closest to failure. This is even more conservative in case of stationary signals: Figure 2 in the main body of the paper provides a heuristic explanation of this and the implemented modification for its improvement, which leads to the addition of step 6 to the procedure.

-Step 6: stationary signal identification and RT(t) estimation as MTTF(t).

This is the key modification with respect to the previously presented procedure [Zio et al., 2009a]. The details are presented in Section 2 of the main body of the paper.


Roisin Lonergan
email: roisin.lonergan@st-vincents.ie

K Kinsella

P Fitzpatrick

J Brady

B Murray

C Dunne

M Duggan

M Mckenna

M Hutchinson

Multiple sclerosis prevalence in Ireland: relationship to vitamin D status and HLA genotype

Keywords: Multiple Sclerosis prevalence, epidemiology, 25(OH) D, vitamin D, latitudinal variation, genetics, HLA

Introduction

Numerous studies support the hypothesis that multiple sclerosis (MS) is a multifactorial autoimmune disease arising from a complex interplay between genetic and environmental factors. The HLA-DRB1*1501 allele appears to be the predominant risk allele in northern European populations [START_REF] Peltonen | Old suspects found guilty-the first genome profile of multiple sclerosis[END_REF] and the HLA DRB1*1501-DQB1*0602 haplotype confers greatest susceptibility to MS in Ireland [START_REF] Dunne | Human leucocyte antigen class II polymorphism in Irish paients with multiple sclerosis[END_REF]. More recently, linkage to other candidate genes has been demonstrated, including IL2-Ra, IL-7 RA [START_REF] Peltonen | Old suspects found guilty-the first genome profile of multiple sclerosis[END_REF], EVI5 (ecotropic viral integration site 5) [START_REF] Ebers | Environmental factors and multiple sclerosis[END_REF] and IFNAR [START_REF] Gilli | Expression and regulation of IFNα/β receptor in IFNβtreated patients with multiple sclerosis[END_REF].

Many environmental factors have been implicated in the aetiology of MS, and may interact with genetic susceptibility. The role of the environment is supported by the distinctive geographical distribution of MS [START_REF] Acheson | Some comments on the relationship of the distribution of multiple sclerosis to latitude, solar radiation, and other variables[END_REF][START_REF] Kurtzke | Epidemiology of multiple sclerosis in U.S. veterans: 1. Race, sex, and geographic distribution[END_REF][START_REF] Miller | Multiple sclerosis in Australia and New Zealand: are the determinants genetic or environmental?[END_REF], with the greatest incidences recorded at high latitudes north and south of the equator. [START_REF] Ascherio | Environmental risk factors for multiple sclerosis, part I: the role of infection[END_REF][START_REF] Pugliatti | The worldwide prevalence of multiple sclerosis[END_REF][START_REF] Munger | Serum 25hydroxyvitamin D levels and risk of multiple sclerosis[END_REF] It is proposed that the latitudinal variation in MS prevalence is inversely related to sunshine exposure and associated vitamin D production. [START_REF] Ebers | Environmental factors and multiple sclerosis[END_REF] The hypothesis regarding the environmental contribution of vitamin D is supported by migration studies, which demonstrate an association between place of early residence and MS. [START_REF] Dean | Age at immigration to England of Asian and Caribbean immigrants and the risk of developing multiple sclerosis[END_REF][START_REF] Hammond | The age-range of risk of developing multiple sclerosis: evidence from a migrant population in Australia[END_REF]Geospatial analysis in North America and continental USA showed a strong negative correlation between MS distribution and ultraviolet B (UV-B) exposure [START_REF] Mckenna | Differences in vitamin D status between countries in young adults and the elderly[END_REF].During winter at high latitudes sun exposure does not generate any vitamin D, thus vitamin D insufficiency is common. [START_REF] Mckenna | Differences in vitamin D status between countries in young adults and the elderly[END_REF] Lower 25(OH)D levels during relapses, in addition to blunted parathyroid hormone (PTH) response [START_REF] Soilu-Hanninen | A longitudinal study of serum 25-hydroxyvitamin D and intact PTH levels indicate the importance pf vitamin D and calcium homeostasis regulation in multiple sclerosis[END_REF], suggest that activated vitamin D may have immunomodulatory effects and influence disease activity [START_REF] Hayes | The immunological function of the vitamin D endocrine system[END_REF]. Activated vitamin D inhibits CD4+T cell and MBP-specific IL-17-secreting T cell proliferation. [START_REF] Correale | Immunomodulatory effects of vitamin D in multiple sclerosis[END_REF] Clinically this is supported by seasonal fluctuation in the number of relapses [START_REF] Australia | Genome-wide association study identifies new multiple sclerosis susceptibility loci on chromosomes 12 and 20[END_REF] and gadolinium-enhancing lesions on MRI. [START_REF] Auer | Seasonal fluctuations of gadolinium-enhancing magnetic resonance imaging lesions in multiple sclerosis[END_REF][START_REF] Embry | Vitamin D and seasonal fluctuations of gadolinium-enhancing magnetic resonance imaging lesions in multiple sclerosis[END_REF] Few studies on the therapeutic effects of vitamin D in RRMS have been performed.

There is no international consensus on optimal 25(OH)D levels or recommended daily intake, and the precise timing of vitamin D action, whether a discrete period or ongoing, remains unknown.

Recent reports suggest that MS prevalence is increasing worldwide. [START_REF] Ebers | Environmental factors and multiple sclerosis[END_REF][START_REF] Alvaro | Temporal trends in the incidence of multiple sclerosis. A systematic review[END_REF][START_REF] Gray | Factors in the rising prevalence of multiple sclerosis in the north-east of Ireland[END_REF] The variation of prevalence with latitude was previously demonstrated in Ireland [START_REF] Mcguigan | Latitudinal variation in the prevalence of multiple sclerosis in Ireland, an effect of genetic diversity[END_REF]which at 53° North lies in a high-risk MS zone where there is minimal skin production of vitamin D between late October and late March due to the attenuation of UV-B rays. Based on previous Irish studies, we know that genetic factors contribute to MS susceptibility [START_REF] Dunne | Human leucocyte antigen class II polymorphism in Irish paients with multiple sclerosis[END_REF], but vitamin D may influence risk. We hypothesized that: 1) the latitudinal gradient in MS prevalence persists in Ireland and 2) greater MS prevalence at higher latitudes relates to lower 25(OH)D levels compared to southern areas and 3) serum 25(OH)D is lower in patients with MS than in healthy controls.

Our aims were to 1) compare the prevalence of MS in three areas in Ireland (two rural areas: counties Donegal (northern latitude) and Wexford (southern latitude), and one urban area in southeast Dublin city (intermediate latitude), 2) determine differences in

Patients and methods

Patients of Irish nationality with MS in each area (Donegal, latitude 54°8'-55°43', Wexford latitude 52°20'-52°44', and South Dublin, latitude 53°17') were identified from multiple sources. The ascertainment sources included neurologists and MS nurse specialists, general practitioners, hospital physicians, Public Health Nurses, the MS Society, pharmacies, Cheshire homes and respite facilities. Control participants, frequency-matched for age and sex, were recruited in each of the three areas. The prevalence date was established as December 31 st , 2007. Informed consent was obtained from each participant. For patients with MS, their MS subtype was recorded [START_REF] Lublin | Defining the clinical course of multiple sclerosis: results of an international survey[END_REF] and the expanded disability status scale (EDSS) was calculated. [START_REF] Kurtzke | Rating neurologic impairment in multiple sclerosis: an expanded disability status scale (EDSS)[END_REF] The diagnosis of MS was based on the revised McDonald criteria [START_REF] Polman | Diagnostic Criteria for Multiple Sclerosis: 2005 Revisions to the "McDonald Criteria[END_REF], and, where necessary, was confirmed by referring to patients' medical notes. Control participants were recruited in each of the three areas and included spouses and relatives of people coefficients of variations (CVs) are 6.2% and 7.7% at concentrations of 28.8 nmol/L and 105.4 nmol/L. The intra-assay CVs are 3.0% and 2.7% at concentrations of 28.9 nmol/L and 73.9 nmol/L. Serum 25(OH)D levels of <50 nmol/L were considered to be vitamin D insufficient, and levels below 25 nmol/L were considered to be deficient. [START_REF] Mckenna | Secondary hyperparathyroidism in the elderly: means to defining hypovitaminosis D[END_REF][START_REF] Holick | Vitamin D deficiency[END_REF] Serum intact PTH was measured by electrochemiluminescent immunoassay on the Elecsys 1010 platform (Roche Diagnostics, Basel, Switzerland).

Intra-assay CV was <5.8% and inter-assay CV was <7.1%.

HLA typing

Genomic DNA was isolated from EDTA-anticoagulated peripheral blood samples using the GenoM-6 (Genovision, Vienna, Austria) DNA isolation and purification system. Low resolution HLA-DRB1/3/4/5 typing was performed using the polymerase chain reaction sequence-specific oligonucleotide probe assay (Invitrogen-Dynal Biotech, Bromborough, UK). HLA types were determined using Dynal RELI-SCAN and RELI SSO PMP software.

Statistical analysis

SAS version 9 (SAS Institute Inc., Cary, NC, USA) was used. The Z test compared rates of MS in two areas, and the Chi square test was used for comparison of rates in three areas. Win PEPI [START_REF] Abramson | WINPEPI (PEPI-for-Windows): computer programs for epidemiologists[END_REF] was used for comparison of rates. For comparison of means and proportions, the student's t test, analysis of variance and chi square test were employed, and Pearson's correlation coefficient was calculated for correlation.

Age standardized rates were calculated using direct standardization, with the national population from the most recent census as the standard population.[27] For all estimates, 95% confidence intervals were calculated.

Results

Prevalence

In total, 632 MS patients were identified in the three regions: 329 in Donegal, 173 in Wexford and 130 in South Dublin (Table 1 2). We could not compare age-standardized prevalence rates over time because this data was not included in the 2001 prevalence study.

Overall, the female to male ratio in MS patients was 2.0:1. In Wexford, the ratio was There was no significant difference between mean 25(OH)D levels in patients in Donegal and Wexford (Donegal 36.9 nmol/L and Wexford 39.7 nmol/L, p = 0.284).

The proportion of MS patients from all three study areas with vitamin D deficiency (<25 nmol/L) (28.3%) was significantly higher than the proportion of control subjects (19.2%) (p = 0.004).

There was a significant inverse correlation between 25(OH)D and PTH levels in MS patients (r= -0.17; p < 0.002) and in controls (r= -0.24; p <0.0003). Vitamin D deficiency (<25 nmol/L) was observed in 93 patients (28.3%) with a mean of 20.

Discussion

In keeping with recent epidemiological studies of MS [START_REF] Ebers | Environmental factors and multiple sclerosis[END_REF][START_REF] Alvaro | Temporal trends in the incidence of multiple sclerosis. A systematic review[END_REF][START_REF] Gray | Factors in the rising prevalence of multiple sclerosis in the north-east of Ireland[END_REF][START_REF] Mcguigan | Latitudinal variation in the prevalence of multiple sclerosis in Ireland, an effect of genetic diversity[END_REF]we have observed an increase in the prevalence of MS nationally between 2001 and 2007, with a significant increase in the most northerly region of Donegal. We confirmed the latitudinal gradient found in 2001 [START_REF] Mcguigan | Latitudinal variation in the prevalence of multiple sclerosis in Ireland, an effect of genetic diversity[END_REF], with a significantly greater MS prevalence in the Northwest compared with the Southeast. This latitudinal gradient has been attributed to genetic factors, based on ancestral differences in each area (Celtic in the Northwest and Anglo-Norman in the Southeast), and the heterogeneity in distribution of HLA class II associations in (specifically the HLA DRB1*1501-DQB1*0602 haplotype) between these areas. The ethnic make-up of the study cohort was relatively homogenous. Although Dublin county has experienced significant migration over the past decade, the particular area studied is not an area of marked migration, and population ethnicity was comparable with both Donegal and Wexford. Increased prevalence may be influenced by numerous factors in addition to increased incidence, but our case ascertainment and diagnostic techniques have not changed greatly since the last prevalence study in 2001. [START_REF] Mcguigan | Latitudinal variation in the prevalence of multiple sclerosis in Ireland, an effect of genetic diversity[END_REF]Incidence could not be calculated due to inability to ascertain year of (MS) onset data in all identified cases. These particular cases were identified from at least two sources, but exact timing of MS onset was not clearly recorded, and patients declined interview.

Our study hypothesis is based on the idea that serum vitamin D levels in adulthood, after MS has presented, are relevant, or at least reflect what has happened previously, a limitation imposed by the nature of our study. Evidence from monozygotic (MZ) twin studies in MS [START_REF] Orton | Evidence for genetic regulation of vitamin D status in twins with multiple sclerosis[END_REF] supports the hypothesis that MS is strongly influenced by environmental factors in genetically susceptible individuals, that vitamin D levels are genetically regulated, and that MS itself is not the cause of low vitamin D (e.g. limited sun exposure imposed by indoor lifestyle). Adoptee [START_REF] Ebers | A genetic basis for familial aggregation in multiple sclerosis. Canadian Collaborative Study Group[END_REF] step-sibling [START_REF] Dyment | Multiple sclerosis in stepsiblings: recurrence risk and ascertainment[END_REF] and half-sibling studies [START_REF] Ebers | Parent-of-origin effect in multiple sclerosis: observations in halfsiblings[END_REF] indicate that the environmental component is ubiquitous, rather than limited to the shared family environment. However, the exact timing of vitamin D action, and whether adult levels reflect those in the past, is as yet unknown.

Chronic low concentrations, seasonal fluctuation, or acute effects may all increase risk, and this limits our ability to draw clear conclusions. Insufficient pre-and perinatal vitamin D levels may imprint on functional characteristics of various tissues, increasing risk of adult onset MS, by interacting with other genetic susceptibility factors [START_REF] Ramagopalan | Expression of the Multiple Sclerosis-Associated MHC Class II Allele HLA-DRB1*1501 Is Regulated by Vitamin D[END_REF] and adult exposures (including hypovitaminosis D during adult life).

MS displays a month-of-birth effect, with an excess of patients being born in spring and a deficit in winter; the effect is most pronounced in familial cases [START_REF] Willer | Timing of birth and risk of multiple sclerosis: population based study[END_REF], implicating interactions between genes and environment related to climate. A recent Scottish study [START_REF] Hannah | Timing of Birth and Risk of Multiple Sclerosis in the Scottish Population[END_REF] showed a marked difference (50%) in MS risk between being born in April versus November, and this month-of-birth association implicates interaction of a seasonal risk factor with loci at or near HLA-DRB1, the region exerting the single strongest genetic effect in MS, during gestation or shortly after birth. [START_REF] Ramagopalan | HLA-DRB1 and month of birth in multiple sclerosis[END_REF] Our study was not sufficiently powered to address month of birth. The nature of this timing of birth effect remains unknown, but maternal vitamin D deficiency is one plausible candidate. [START_REF] Willer | Timing of birth and risk of multiple sclerosis: population based study[END_REF][START_REF] Looker | Serum 25hydroxyvitamin D status of adolescents and adults in two seasonal subpopulations from NHANES III[END_REF][START_REF] Ramagopalan | Expression of the Multiple Sclerosis-Associated MHC Class II Allele HLA-DRB1*1501 Is Regulated by Vitamin D[END_REF] Although monozygotic twins (MZ) share a common genotype, epigenetic differences may cause phenotypic discordance (including disease susceptibility). Examination of differences in DNA methylation and histone acetylation in a large cohort of monozygotic twins [START_REF] Fraga | Epigenetic differences arise during the lifetime of monozygotic twins[END_REF] showed that, although epigenetically indistinguishable during early years of life, marked differences in content and genomic distribution of 5methylcytosine DNA and histone acetylation in older twins affected their geneexpression portrait. The divergence of epigenetic modification patterns with age may be influenced by intrinsic and extrinsic factors such as smoking and diet, and epigenetic drift may occur with aging itself. Fraga et al showed that approximately one third of twins studied displayed epigenetic differences in DNA methylation and histone modification, distributed throughout their genomes and impacting on gene expression. These markers were more distinct in older twins with different lifestyles, highlighting the role of environmental factors in modifying a common genotype to a different phenotype (e.g. discordant MS onset in MZ twins). Environmental modulation of genetic information, by altering the pattern of epigenetic modifications, supports the environment (vitamin D)-genetic interaction findings in our study, and suggests that adult exposures, including vitamin D levels, may play a role in both MS risk and disease modification, even if the critical exposure period began prenatally.

Despite our hypothesis that the prevalence gradient would relate to lower 25(OH)D in MS patients, 25(OH)D levels between the north (Donegal) and south (Wexford) were not significantly different and did not differ from control subjects: winter levels were low in all participants. It is important to remember that optimal levels remain undefined (and may be much higher than current laboratory reference ranges, e.g. ≥100 nmol/L) [START_REF] Bischoff-Ferrari | Estimation of optimal serum concentrations of 25-hydroxyvitamin D for multiple health outcomes[END_REF]. Our view of sufficiency is arbitrary, and we cannot assume that all individuals require the same levels. However, mean winter levels in our Irish cohort were low overall (<50 nmol/L) and significantly more MS patients than controls had very low levels (<25 nmol/L). As previously described by others [START_REF] Soilu-Hanninen | A longitudinal study of serum 25-hydroxyvitamin D and intact PTH levels indicate the importance pf vitamin D and calcium homeostasis regulation in multiple sclerosis[END_REF], there was a blunted PTH response only in MS patients with vitamin D deficiency (<25nmol/L), but its meaning is unclear. Only a small percentage took >1,000IU vitamin D3 daily; maximum intake was <1,500IU, but it is likely that at least 4,000IU daily are required to achieve levels of 100 nmol/L [START_REF] Vieth | Vitamin D supplementation, 25-hydroxyvitamin D concentrations, and safety[END_REF]. Despite differences in solar radiation (UV-B) exposure between areas; the northwest of Ireland experiences on average 319,875 joules/cm 2 solar radiation annually, compared with 331,683 joules/cm 2 solar radiation in the southeast [44]; it is probable that the latitudinal range in a small country such as Ireland is insufficient to demonstrate such an effect. Higher 25(OH)D in southeast Dublin may be genetically regulated, or relate to more widespread use of vitamin D supplements (39%) compared with Wexford (24%), and greater availability and consumption of fortified milk in Dublin. [START_REF] Mckenna | Safety and efficacy of increasing wintertime vitamin D and calcium intake by milk fortification[END_REF] Differences in socio-economic group distribution between these areas may influence supplementation practices. [27] However, in 52% of patients and 55% of controls taking vitamin D3 supplements the dose of ≤ 400IU/day was unlikely to alter serum levels significantly. In addition to the previously discussed limitation imposed by the nature of our study (i.e. hypothesis based on the idea that serum vitamin D levels in adulthood are relevant in MS), we acknowledge that inclusion of spouses as controls is not always ideal, as they are usually in the same environment as affected cases, and may modify their risk by taking vitamin D supplements more readily. However, we found that significantly more patients (18.2%) than controls (11.5%) were taking vitamin D supplements (≥800IU) (p =0.04), despite this potential spousal influence. Spousal controls were included because they were more willing to give samples than people without affected family or partners.

serum 25 (

 25 OH)D and PTH levels between: (a) MS patients in the three regions and (b) MS patients and healthy control subjects, and 3) examine the frequency of HLA DRB1*15 in people with MS in northern and southern counties.

  with MS in addition to unrelated volunteers. Data from the most recent National Census [27] were used to calculate prevalence. Blood samples from patients and controls were taken in winter months from December to March for serum 25(OH)D and serum PTH measurement, and for DNA extraction. A questionnaire regarding environmental factors including past and current vitamin D supplementation and sun exposure was completed.[28] The Ethics Committee at St. Vincent's University Hospital approved the study. Serum 25(OH)D and serum PTH analysis Serum 25(OH)D was measured by a competitive radioimmunoassay (Immunodiagnostic Systems Limited, Boldon, Tyne & Wear, UK). Inter-assay

  Age standardised MS prevalence was: Donegal -all MS patients 257.8/10 5 (95% CI 253.2-262.4), male MS patients 157.6/10 5 (95% CI 152.7-162.6), female MS patients 358.3/10 5 (95% CI 350.3-366.5); Wexford-all MS patients 152.0/10 5 (95% CI 148.4-155.6), male MS patients 119.0/10 5 (95% CI 114.8-123.4), female MS patients 176.0/10 5 (95% CI 170.4-181.7); southeast Dublin-all MS patients 123.2/10 5 (95% CI 120.0-126.4), male MS patients 84.2/10 5 (95% CI 80.6-87.9), female MS patients 160.6/10 5 (95% CI 155.2-166.1). The 3 counties were significantly different from each other when comparing age standardized prevalence in all MS patients (p <0.001), males with MS (p <0.001) and females with MS (p <0.001). In Donegal, MS prevalence had increased significantly since 2001; prevalence was 184.6/10 5 in 2001, (95%CI 162.0-209.5) and 290.3/10 5 in 2007, (95%CI 238.7-255.5) (p=0.001). There was a non-significant increase in prevalence in Wexford; prevalence was 120.7/10 5 in 2001, (95%CI 100.6 to 143.8) and 144.8/10 5 in 2007, (95%CI 138.3 -151.9) (p= 0.117). Prevalence was lowest in southeast Dublin (130.0/10 5 in 2007) but we did not have previous data for comparison (Table

  3 nmol/L (range 13.3-24.6) and a mean PTH of 44.8 ng/L (range 16.7-86.3). There was a non-significant inverse correlation between 25(OH)D and PTH in this group of patients (p =0.28). Vitamin D deficiency (<25 nmol/L) was observed in 44 control subjects (19.2%) with a mean of 20.8 nmol/L (range 14.3-24.6 nmol/L) and a mean PTH of 48.3ng/L (range 5.0-131.5 ng/L). A significant inverse correlation between 25(OH)D and PTH in the control subjects with 25(OH)D <25 nmol/L (r= -0.47; p =0.0013) indicated a normal response to vitamin D deficiency. Overall, significantly more patients (18.2%) than controls (11.5%) were taking vitamin D supplementation of >800IU/day (800IU-1500IU/day) (p =0.04), but there was no significant difference in proportions of patients (5.7%) and controls (8.6%) taking vitamin D supplementation of 1,300IU-1,500IU/day (p >0.05). None of the cohort took >1,500IU vitamin D supplementation/day. We stratified serum 25(OH)D results, looking separately at patients and controls not taking vitamin D supplementation. Mean serum vitamin D level was similar in patients (38.4 nmol/L) and controls (37.7 nmol/L) after stratification (p >0.05). There was no significant difference in mean 25(OH)D levels between patients (41.0 nmol/L) and controls (38.4 nmol/L) who were taking up to 800IU vitamin D/day (p >0.05). The difference was significant only in those with very low levels ('deficiency'= ≤25 nmol/L), as before, with a greater proportion of MS patients having very low levels (p = 0.004). Mean serum 25(OH)D levels overall (patients and controls) were significantly lower in those taking no supplements or <800IU vitamin D/day (38.1 nmol/L) than in those taking >800IU/day (53.1 nmol/L) (p < 0.001). Mean serum 25(OH)D were also significantly lower in those taking no supplements or <1000IU vitamin D/day (39.9 nmol/L) than in those taking 1,000-1,500IU/day (47.9 nmol/L) (p = 0.03). Vitamin D supplement data are displayed in table 4. Although we did not observe significant relationships between serum 25(OH)D levels and age at MS onset (p =0.26) or the presence of relapse in past 2 years (p =0.7) , our data were underpowered to address questions regarding these factors, or to comment on any association between MS risk and season of birth, early sun exposure (two or more hours of daily summer sun exposure in childhood (p= 0.21) or vitamin D supplementation (current or up to age 15, p= 0.76). Frequency of HLADRB1*15 allele The HLADRB1*15 allele was present in 62.9% of MS patients, compared with 35.2% of control subjects (p < 0.01). The frequency was significantly higher in Donegal patients (69%) compared with South Dublin patients (52%) (p < 0.04). Homozygosity for the HLADRB1*1501-DQB1*0602 haplotype occurred significantly more frequently in Donegal compared with Wexford (p = 0.032).

  The HLA DRB1*15 allele, associated with MS susceptibility, was significantly overrepresented in MS patients compared with controls, and HLA-DRB1*15 frequency was significantly lower in southeast Dublin compared with areas of higher MS prevalence. Active vitamin D (1, 25-dihydroxyvitamin D) binds the intranuclear vitamin D receptor (VDR), whose gene locates to chromosome 12q13.1. The resulting complex binds the retinoid X receptor (RXR), resulting in a VDR-RXR complex, which then binds the vitamin D response element (VDRE).[START_REF] Smolders | Vitamin D as an immune modulator in multiple sclerosis, a review[END_REF] Recently, a VDRE was found on HLA DRB1*1501, the major MS susceptibility allele.[47] Insufficient HLA DRB1*1501 expression -in the setting of low vitamin D and suboptimal binding to the VDRE, with poor up regulation of HLA DRB1*1501 -could affect central deletion of auto-reactive T cells. Thus people possessing HLA DRB1*1501 may be at high risk of auto-immunity (MS) in the setting of vitamin D deficiency and this combination may explain the high prevalence of MS in the Irish population. The latitudinal variation in MS in Ireland is not explained by 25(OH)D levels alone, but most likely the interaction with genetic susceptibility. In a gene-environment interaction, lacking either component can prevent clinical expression of disease, leaving one factor apparently unrelated to risk when both are required. Therefore, similar overall serum 25(OH)D levels in patients and controls does not outrule the role of low vitamin D in MS pathogenesis in our cohort.

  In conclusion, although unproven, one might hypothesise that combined high frequency of HLADRB1*15 and widespread vitamin D insufficiency may explain the high, increasing MS prevalence in the Irish population. Vitamin D may only be one factor in a common mechanism by which environment interacts, but it is a potentially modifiable risk. Vitamin D insufficiency is common throughout the population, an important public health observation. The MS prevalence gradient in Ireland is not explained by variation in 25(OH)D levels alone, but differences in frequency of the MS risk allele HLADRB1*1501 and risk haplotype HLADRB1*1501-DQB1*0602 combined with low serum vitamin D levels probably influence latitudinal variation.

Figure 1

 1 Figure 1 Map of Ireland demonstrating the 3 areas studied and their latitudes: Donegal

  ). MS prevalence on December 31 st , 2007 was: Donegal 290.3/10 5 , (95% CI; 238.7 -255.5), Wexford 144.8/10 5 , (95% CI; 138.3 -151.9), southeast Dublin 127.8/10 5 , (95% CI; 131.6 -146.7). Prevalence was significantly higher in Donegal than Wexford (p <0.0001). Although prevalence in Wexford was higher than in South Dublin this did not reach significance (p =0.72).

  OH)D and PTH levels were determined in 329 MS patients (61.3% of all ascertained patients), and 226 healthy control subjects (Table3).

	25(OH)D and PTH levels
	25(The mean level of serum 25(OH) D was similar in MS patients (38.6 nmol/L, SD
	22.5) and controls (36.4 nmol/L, SD 16.4 (p >0.05)). Mean 25(OH)D levels in
	patients in Donegal were 36.9 nmol/L (range: 13-97 nmol/L), in Wexford 39.7
	nmol/L (15-114 nmol/L) and in southeast Dublin 50.7 nmol/L (13-161 nmol/L). Mean
	levels were significantly higher in southeast Dublin, the area with lowest MS
	prevalence, than in Donegal (p <0.0001) and Wexford (p <0.0001).
	1.5:1 (2.2 in 2001), in Donegal 2.1:1 (2.48 in 2001) and in southeast Dublin 2.6:1 (no
	earlier data).

Table 1 : Demographic data of MS patients ascertained in each region of Ireland on 31/12/2007

 1

		Donegal	Wexford South Dublin	Overall
	Mean age (years)	46.7	47.9	49.5	48.0
	Range	(21-80)	(20-75)	(27-89)	(20-89)
	F: M ratio	2.1:1	1.5:1	2.6:1	2.0:1
	MS duration (years)	15.3	14.9	17.2	15.8
	Range	(1-54)	(2-49)	(1-53)	(1-54)
	Relapsing remitting %	53.3	47.0	53.1	51.1
	Secondary progressive %	37.6	43.2	39.2	39.7
	Primary progressive %	9.1	9.8	7.7	9.2

Mean age MS onset (yrs) Range 31.6 (8-64) 32.5 (14-54) 32.8 (16-61) 32.3 (8-64) Mean EDSS Range

Table 2 : Population and MS prevalence data in Donegal, Wexford and South Dublin in 2001 and 2007 Region (Latitude)

 2

						Age	Age	Age
		Census 1996*	MS cases 2001	Census 2006*	MS cases 2007	standardized prevalence rates:	standardized prevalence rates:	standardized prevalence rates:
						All MS cases	Male	Female
			240		329
	Donegal (54° 8'-	129,994	Prevalence: 184.6/10 5	113,347	Prevalence: 290.3/10 5
	55°43')		(95% CI:		(95% CI:
			162.0-209.5)		238.7

-255.5) p <0.001 257.8/10 5 253.2-262.4 157.6/10 5 152.7-162.6 358.3/100,000 350.3-366.5 Wexford

			126
			Prevalence:
	(52° 2'-	104,372	120.7/10 5	119,442
	52°44')		(95% CI:
			100.6-143.8)

Table 4 : Vitamin D supplement doses taken by MS patients and controls

 4

		38.57	36.41	36.90	39.71	50.74
	Range	13.30-161.90	14.30-97.20 13-97	15-114	13-161
	Mean PTH (ng/L)	41.07	42.04	42.93	39.94	38.95
	Range	14.90-321.70	5-131.50	13-97	18-90	8-39
	% with HLA-DRB1*1501	62.90%	35.20% (P<0.01)	69%	63%	52%

Legend-IU international units; % percentage Of note, standard multivitamins, usually containing only 400 IU/day, are unlikely to alter serum levels of vitamin D. Legend for figure 1

Acknowledgements

The authors would like to sincerely thank the consultant neurologists at AMNCH, Tallaght, St. James' Hospital, Beaumont Hospital and Mater Misericordiae hospitals in Dublin, Sligo General Hospital and the Western Regional Hospital, Limerick and Waterford Regional Hospital, and also the general practitioners in each region, for granting access to data of MS patients in our study areas, and also Dr. Bannan, Consultant Physician at Letterkenny Hospital and Dr. Quigley Consultant Physician at Wexford General Hospital, who allowed us to use their clinic rooms for weekend research clinics. We are also grateful to MS Ireland for much help during the study, and for ongoing support. We also thank the patients with MS and the control volunteers who agreed to participate in this study.

Licence for Publication:

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd to permit this article (if accepted) to be published in JNNP and any other BMJPGL products and sublicences such use and exploit all subsidiary rights, as set out in our licence.

Competing Interest: None declared.