This paper proposes a clustered model reduction method for semistable positive linear systems evolving over directed networks. In this method, we construct a set of clusters, i.e., disjoint sets of state variables, based on a notion of cluster reducibility, defined as the uncontrollability of local states. By aggregating the reducible clusters with aggregation coefficients associated with the Frobenius eigenvector, we obtain an approximate model that preserves not only a network structure among clusters, but also several fundamental properties, such as semistability, positivity, and steady state characteristics. Furthermore, it is found that the cluster reducibility can be characterized for semistable systems based on a projected controllability Gramian that leads to an a priori H-2-error bound of the state discrepancy caused by aggregation. The efficiency of the proposed method is demonstrated through an illustrative example of enzyme-catalyzed reaction systems described by a chemical master equation. This captures the time evolution of chemical reaction systems in terms of a set of ordinary differential equations. (C) 2015 Elsevier Ltd. All rights reserved.