function [N,E] = rentian_scaling_3d(A,XYZ,n,tol) % RENTIAN_SCALING_3D Rentian scaling for networks embedded in three dimensions. % % [N,E] = rentian_scaling_3d(A,XYZ,n,tol) % % Physical Rentian scaling (or more simply Rentian scaling) is a property % of systems that are cost-efficiently embedded into physical space. It is % what is called a "topo-physical" property because it combines information % regarding the topological organization of the graph with information % about the physical placement of connections. Rentian scaling is present % in very large scale integrated circuits, the C. elegans neuronal network, % and morphometric and diffusion-based graphs of human anatomical networks. % Rentian scaling is determined by partitioning the system into cubes, % counting the number of nodes inside of each cube (N), and the number of % edges traversing the boundary of each cube (E). If the system displays % Rentian scaling, these two variables N and E will scale with one another % in loglog space. The Rent's exponent is given by the slope of log10(E) % vs. log10(N), and can be reported alone or can be compared to the % theoretical minimum Rent's exponent to determine how cost efficiently the % network has been embedded into physical space. Note: if a system displays % Rentian scaling, it does not automatically mean that the system is % cost-efficiently embedded (although it does suggest that). Validation % occurs when comparing to the theoretical minimum Rent's exponent for that % system. % % INPUTS: % % A: MxM adjacency matrix. % Must be unweighted, binary, and symmetric. % XYZ: Matrix of node placement coordinates. % Must be in the form of an Mx3 matrix [x y z], where M is % the number of nodes and x, y, z are column vectors of node % coordinates. % n: Number of partitions to compute. Each partition is a data % point. You want a large enough number to adequately % estimate the Rent's exponent. % tol: This should be a small value (for example 1e-6). % In order to mitigate the effects of boundary conditions due % to the finite size of the network, we only allow partitions % that are contained within the boundary of the network. This % is achieved by first computing the volume of the convex % hull of the node coordinates (V). We then ensure that the % volume of the convex hull computed on the original node % coordinates plus the coordinates of the randomly generated % partition (Vnew) is within a given tolerance of the % original (i.e. check abs(V - Vnew) < tol). Thus tol, should % be a small value in order to make sure the partitions are % contained largely within the boundary of the network, and % thus the number of nodes and edges within the box are not % skewed by finite size effects. % % OUTPUTS: % % N: nx1 vector of the number of nodes in each of the n partitions. % E: nx1 vector of the number of edges crossing the boundary of % each partition. % % Subsequent Analysis: % % Rentian scaling plots are created by: figure; loglog(E,N,'*'); % % To determine the Rent's exponent, p, we need to determine % the slope of E vs. N in loglog space, which is the Rent's % exponent. There are many ways of doing this with more or less % statistical rigor. Robustfit in MATLAB is one such option: % % [b,stats] = robustfit(log10(N),log10(E)) % % Then the Rent's exponent is b(1,2) and the standard error of the % estimation is given by stats.se(1,2). % % Note: n=5000 was used in Bassett et al. 2010 in PLoS CB. % % Reference: % Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg, % Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient % physical embedding of topologically complex information processing % networks in brains and computer circuits. PLoS Comput Biol, 2010, % 6(4):e1000748. % % Modification History: % % 2010: Original (Dani Bassett) % Dec 2016: Updated code so that both partition centers and partition % sizes are chosen at random. Also added in a constraint on % partition placement that prevents boxes from being located % outside the edges of the network. This helps prevent skewed % results due to boundary effects arising from the finite size % of the network. (Lia Papadopoulos) % % determine the number of nodes in the system M = numel(XYZ(:,1)); % rescale coordinates so that they are all greater than unity XYZn = XYZ - repmat(min(XYZ)-1,M,1); % compute the area of convex hull (i.e. are of the boundary) of the network [~,V] = convhull(XYZn(:,1),XYZn(:,2),XYZn(:,3)); % min and max network coordinates xmin = min(XYZn(:,1)); xmax = max(XYZn(:,1)); ymin = min(XYZn(:,2)); ymax = max(XYZn(:,2)); zmin = min(XYZn(:,3)); zmax = max(XYZn(:,3)); % initialize vectors of number of nodes in box and number of edges crossing % box N = zeros(n,1); E = zeros(n,1); % create partitions, and count the number of nodes inside the partition (N) % and the number of edges traversing the boundary of the partition (E) nPartitions = 0; % create partitions, and count the number of nodes inside the partition (N) % and the number of edges traversing the boundary of the partition (E) while nPartitions<(n+1) % variable to check if partition center is within network boundary % OK if inside == 1 inside = 0; while inside == 0 % pick a random (x,y,z) coordinate to be the center of the box randx = xmin+(xmax-xmin)*rand(1); randy = ymin+(ymax-ymin)*rand(1); randz = zmin+(zmax-zmin)*rand(1); % make sure the point is inside the convex hull of the network newCoords = [XYZn; [randx randy randz]]; [~,Vnew] = convhull(newCoords(:,1),newCoords(:,2),newCoords(:,3)); % if the old convex hull area and new convex hull area are equal % then the box center must be inside the network boundary. if isequal(V,Vnew)==0 inside = 0; else inside = 1; end end % determine the approximate maximum distance the box can extend, given % the center point and the bounds of the network deltaX = min(abs(xmax-randx),abs(xmin-randx)); deltaY = min(abs(ymax-randy),abs(ymin-randy)); deltaZ = min(abs(zmax-randz),abs(zmin-randz)); deltaLmin = min([deltaX deltaY deltaZ]); % variable to check if partition is within network boundary % OK if inside == 1 inside = 0; while inside == 0 % pick a random (side length)/2 that is between 0 and the % max possible deltaL = deltaLmin*rand(1); % (x,y,z) coordinates for corners of box boxCoords = [randx - deltaL randy - deltaL randz - deltaL; ... randx - deltaL randy - deltaL randz + deltaL; ... randx - deltaL randy + deltaL randz - deltaL; ... randx - deltaL randy + deltaL randz + deltaL; ... randx + deltaL randy - deltaL randz - deltaL;... randx + deltaL randy - deltaL randz + deltaL; ... randx + deltaL randy + deltaL randz - deltaL; ... randx + deltaL randy + deltaL randz + deltaL]; % check if all corners of box are inside the convex hull of the % network newCoords = [XYZn; boxCoords]; [~,Vnew] = convhull(newCoords(:,1),newCoords(:,2),newCoords(:,3)); % make sure the new convex hull that includes the partition corners % is within a certain tolerance of the original convex hull area. if abs(V-Vnew)>tol inside = 0; else inside = 1; end end % Find nodes inside the box, edges crossing the boundary L = find(XYZn(:,1)>(randx-deltaL) & XYZn(:,1)<(randx+deltaL) ... & XYZn(:,2)>(randy-deltaL) & XYZn(:,2)<(randy+deltaL) ... & XYZn(:,3)>(randz-deltaL) & XYZn(:,3)<(randz+deltaL)); if ~isempty(L) == 1 nPartitions = nPartitions+1; % count edges crossing the boundary of the cube E(nPartitions,1) = sum(sum(A(L,setdiff(1:M,L)))); % count nodes inside of the cube N(nPartitions,1) = numel(L); end end return;